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COMPUTATION  OF  LMS  AND  MATCHED  DIGITAL  FILTERS 
FOR  OPTICAL  CLUTTER  SUPPRESSION 


I.  INTRODUCTION 

Matched  spatial  filters  are  commonly  used  for  suppressing  clutter  and  detecting  signals  in  either 
real-time  or  imaged  output  of  electro-optic  sensors.  Least-mean-square  (LMS)  spatial  filters  have 
been  so-used  in  infrared  (IR)  applications.  This  article  extends  certain  known  methods  of  calculating 
digital  impulse  responses  for  these  filters  so  that  the  methods  are  applicable  with  both  one-  and  two- 
dimensional  (1-D  and  2-D)  forms  of  rather  general  signal  and  clutter  models.  Concomitantly,  2-D 
LMS  filters  are  derived  for  the  first  time,  and  their  make-up  and  operation  are  investigated  and 
explained.  A  primary  objective  of  this  article  is  to  enable  readers  to  compute  digital  impulse-response 
weights  of  an  LMS  or  matched  filter  for  any  signal  and  clutter  representable  by  models  of  the 
assumed  types.  Accordingly,  the  design  methods  are  described  fully  and  illustrated  with  examples  so 
that  they  can  be  applied  easily  in  diverse  conditions. 

The  methods  can  be  used  with  signals  from  scanning  or  staring  sensors  viewing  extended  or 
point  sources  against  variable  backgrounds,  but  signal  shape  and  orientation  must  always  be  known. 
Since  both  point  and  extended  sources  have  finite  images,  a  signal  from  either  can  be  resolved  (con¬ 
sist  of  more  than  one  nonzero  sample).  The  signal  shape  then  varies  depending  on  where  the  image 
&lls  in  relation  to  the  sample  sites.  In  addition,  an  extended  source  eclipses  the  background,  which 
results  in  variable  source  contrast  and  further  signal  variability.  Thus  the  design  methods  are  useful  if 
signal  shape  is  representable  by  an  average  or  if  detection  of  a  particular  signal  is  imperative.  If  the 
signal  is  extremely  variable  and  must  be  detected  in  ail  its  guises,  the  present  methods  fail.  For  rea¬ 
sons  to  be  explained,  the  matched-filter  design  technique  requires  knowledge  of  cluner  parameters, 
whereas  the  LMS  design  technique  does  not.  though  qualitative  knowledge  of  the  background  may  be 
useful  for  choosing  the  LMS  clutter  model.  Needed  information  about  signal  and  clutter  is  assumed 
to  be  available  from  images  of  sensor  output,  obtained  either  experimentally  or  theoretically. 

Matched  filters  are  far  older  and  better  known  than  LMS  filters,  having  been  invented  by  D.  O. 
North  in  1943.'  The  first  matched  filter  was  designed  for  detection  of  radar  signals  in  white  noise. 
Subsequently  the  matched-filter  concept  was  extended  to  other  types  of  noise  and  to  optical  sig¬ 
nals.2-5  The  first  matched  filters  were  for  continuous  input  and  were  difficult  to  make.  Develop¬ 
ment  of  surface-acoustic-wave  delay  lines  and  digital  electronics  overcame  the  manufacturing  prob¬ 
lems  and  greatly  increased  the  use  of  matched  filters  for  processing  both  continuous  and  sampled 
input.6  Today.  2-D  digital  matched  filters  are  widelv  used  for  pattern  recognition  and  image  process¬ 
ing.78 

Matched  filters  result  from  maximizing  the  ratio  of  filtered  signal  energy  to  average  filtered 
noise  power.  There  are  various  ways  to  derive  digital  filter  weights  from  this  criterion.  The  method 
elaborated  here  was  introduced  for  IR  signals  by  Otazo.  Parenti.  and  Tung,  who  employed  an  aniso¬ 
tropic  clutter  model  based  on  analysis  of  measured  IR  backgrounds  '4  10  Their  work  resulted  in  the 
now  well  known  fourth-derivative  filters,  obtained  by  estimating  the  mixed  partial  derivative 
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(5*j  /dx2dy2)  at  sampled  points  of  the  signal  s.  Kay  and  others,  who  also  analyzed  measured  back¬ 
grounds,  have  suggested  that  in  some  respects  an  isotropic  clutter  model  may  be  preferable.11'13  In 
this  article,  the  design  technique  of  Otazo  et  al.  is  extended  by  using  a  generalized  form  of  Kay's  iso¬ 
tropic  model.  The  results  are  presented  as  operator  equations  for  the  filter  weights  and  a  table  of 
operator  arrays. 

Recently  W.  N.  Peters  derived  theoretical  expressions  for  matched-filter  transfer  functions  by 
maximizing  the  signal-to- noise  ratio  (SNR)  of  the  entire  electro-optical  system  from  the  background 
and  source  through  the  digital  filter. 14  The  2-D  filters  obtained  in  this  way  are  called  type  A  or  type  B 
depending  on  whether  signal  and  clutter  are  assumed  to  be  known  in  the  object  plane  or  at  the 
sensor’s  output;  1-D  filters  designed  from  sensor  output  are  called  type  C.  The  matched  filters  of  the 
present  article  are  Peters'  type  B  and  type  C  since  they  are  based  on  sensor  output. 

LMS  filters  are  obtained  from  least-squares  fits  of  signal-plus-clutter  models  to  sensor  output. 
In  fact,  an  LMS  filter  implements  a  regression  analysis  of  the  output  and  thereby  estimates  the 
amount  of  model  signal  present.  Takken  et  al.  introduced  the  1-D  LMS  filter  for  continuous  signals 
and  extended  the  continuous  results  to  sampled  data. 13  Here  analysis  of  the  discrete  case  is  simplified 
by  applying  the  least-squares  method  directly  to  the  sampled  data.  Equations  are  derived  for  the 
weights  of  1-D  LMS  filters  based  on  first-  through  fifth-degree  clutter  polynomials.  The  design  tech¬ 
nique  for  discrete  2-D  LMS  filters  is  then  worked  out.  It  follows  the  same  pattern  as  the  1-D  tech¬ 
nique,  but  there  are  more  equations  with  more  terms.  Equations  are  given  for  the  weights  of  2-D 
LMS  filters  based  on  linear  through  cubic  clutter  models.  Equations  for  higher-order  2-D  filters  are 
easily  derived  but  are  not  given  because  they  are  too  lengthy. 

Frequency-domain  characteristics  of  LMS  filters  receive  considerable  attention  in  reference  IS. 
There  in  Appendix  A  it  is  shown  that  1-D  LMS  filters  are  fixed  bandpass  filters  having  a  special 
property— the  ability  to  block  certain  waveforms  whose  complexity  increases  with  the  degree  of  the 
clutter  polynomial.  Frequency-domain  characteristics  of  LMS  filters  are  not  further  investigated  in 
this  article  because  they  are  not  needed  for  filter  design  and  are  adequately  treated  in  reference  IS. 
Moreover,  emphasizing  the  space-domain  analysis  calls  attention  to  properties  of  LMS  filters  that  are 
not  apparent  in  the  frequency  domain  and  are  not  mentioned  in  reference  IS. 

Because  of  a  similarity  in  names,  it  is  advisable  to  compare  and  contrast  LMS  filters  with  LMS 
adaptive  filters.  The  latter  devices  employ  an  LMS  algorithm  that  adjusts  their  weights  to  minimize 
the  difference  between  the  filter  output  ana  an  externally  supplied  desired  response.16  In  the  present 
context,  the  signal  shape  and  orientation  are  the  desired  response  The  LMS  algorithm  does  not  use 
the  method  of  least  squares  directly  but  rather  approximates  the  solution  of  the  Wiener-Hopf  least- 
squares  equations,  hence  the  names  of  the  algorithm  and  filter.17  The  LMS  algorithm  is  known  to 
converge  with  uncorrelated  and  correlated  stationary  input,  and  with  some  correlated  nonstationary 
input,  but  it  is  not  known  to  converge  unconditionally.1718  Since  the  algorithm  acts  iteratively  and  in 
effect  employs  statistical  samples  of  limited  size,  the  filter  exhibits  settling  time  and  misadjustment 
effects.16-17  These  could  lead  to  problems  with  rapidly  varying  nonstationary  backgrounds.  The  LMS 
filter  discussed  here  is  a  much  simpler  device  whose  weights  are  fixed.  It  nonetheless  conforms  to 
the  input  because  the  parameters  of  its  signal-plus-clutter  model  vary  spontaneously  to  fit  the  back¬ 
ground  and  signal,  if  any,  within  the  filter.  The  quality  of  the  fit  depends  on  the  model's  fidelity  to 
the  local  background  and  so  varies  over  the  scene.  Since  the  LMS  filter  solves  least-squares  regres¬ 
sion  equations  exactly,  there  are  no  questions  of  convergence,  settling  time,  and  misadjustment. 
These  features  make  the  LMS  filter  well  suited  for  suppression  of  optical  clutter,  which  can  ary 
rapidly  in  almost  any  fashion. 

The  best  known  device  derived  by  a  least-squares  analysis  is  the  Wiener  filter.  Its  weights  are 
the  exact  solution  of  the  Wiener-Hopf  equations— the  solution  approximated  by  the  LMS  algorithm. 
Wiener  filters  and  LMS  filters  differ  in  the  information  needed  to  determine  the  filter  weights  and  in 


the  type  of  output.  To  determine  the  weights  of  a  Wiener  filter,  it  is  necessary  to  know  the  auto¬ 
correlation  of  the  input  and  its  cross  correlation  with  the  signal. 19  To  determine  LMS-filter  weights, 
only  signal  and  clutter  models  are  required.  No  statistical  properties  of  the  input  are  needed.  The 
output  of  a  Wiener  filter  estimates  the  signal  itself.  The  output  of  an  LMS  filter  estimates  a  signal 
parameter,  the  amplitude. 

Finally,  a  few  words  about  the  organization  of  this  article  will  be  helpful.  Matched  filters.  1-D 
LMS  filters,  and  2-D  LMS  filters  are  discussed  separately  in  sections  (Q,  m,  IV),  each  having  sub¬ 
sections  A  and  B.  Subsections  A  give  the  basic  analyses  required  for  calculating  filter  weights,  with 
numerical  examples  that  can  be  worked  through  to  test  understanding  of  the  analyses.  Subsections  B 
discuss  the  equations  and  filters  and  analyze  LMS  filters  for  symmetric  signals,  which  present  special 
problems.  Those  who  only  want  to  calculate  filter  weights  can  omit  subsections  B  unless  they  have 
trouble  with  their  computations.  In  that  case  they  may  wish  to  read  these  subsections  and  will  be  in  a 
better  position  to  do  so.  Three  appendices  give  detailed  information  about  matched-filter  clutter 
models,  convolution  arrays,  and  construction  of  these  arrays  to  approximate  derivative  operators. 
Those  who  wish  to  calculate  filter  weights  for  conditions  different  from  the  ones  assumed  in  sections 
H  IV  may  need  to  consult  the  appendices. 

n.  l-D  AND  2-D  MATCHED  FILTERS 

Matched  filters  can  be  derived  for  sampled  or  continuous  sensor  output  by  using  either  clutter 
data  from  a  single  background  or  a  clutter  model  that  approximates  data  from  a  large  set  of  back¬ 
grounds.9,10  The  data-based  and  model-based  approaches  are  complementary.  A  matched  filter  based 
on  data  from  a  single  background  has  the  greatest  signal-to-ciutter  ratio  (SCR)  achievable  with  a 
linear  filter,  but  only  for  the  data  employed.  A  matched  filter  based  on  a  clutter  model  is  more 
widely  applicable;  it  has  nearly  maximum  SCR  against  all  backgrounds  reasonably  well  represented 
by  the  model.  The  performance  of  a  model-based  filter  against  a  specified  background  can  be 
evaluated  by  comparison  with  the  data-based  filter  derived  specifically  for  that  background.  In  this 
way,  other  investigators  have  found  that  matched  filters  based  on  Eq.  (A6)  of  Appendix  A  have  SCRs 
not  more  than  4  dB  below  the  SCRs  of  the  best  possible  matched  filters  against  those  backgrounds.9 
Consequently,  the  matched  filters  developed  in  the  present  work  are  based  on  power-spectral  clutter 
models  somewhat  more  general  than  Eq.  (A6).  Continuous  filters  are  derived  first  and  then  converted 
to  digital  form.  The  impulse  responses  rather  than  the  transfer  functions  are  digitized  because  the 
‘filters  of  interest  are  small— they  act  on  no  more  than  about  50  data  at  once.  In  such  cases  the  filters 
are  more  efficiently  implemented  in  the  space  domain  than  in  the  frequency  domain  (ref.  8.  sec. 
11.3). 


(A)  Analysis  and  Example 

The  one-  and  two-dimensional  filters  will  be  developed  jointly  so  that  the  simpler  1-D  analysis 
will  clarify  the  more  intricate  2-D  treatment.  The  first  step  is  to  derive  the  filter  transfer  functions  by 
maximizing  the  signal-to-average-clutter  ratios.  The  derivations  are  not  given  here  because  they  can 
be  found  in  many  readily  available  sources,  c.g..  refs.  7,8.20-22.  For  present  purposes  the  essential 
results  are  the  equations  of  the  transfer  functions  and  impulse  responses  together  with  definitions  of 
the  quantities  involved.  The  1-D  and  2-D  transfer  functions  are 
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arbitrary  constants; 
imaginary  unit,  j 2  =  -1; 

co-ordinates  of  maximum  signal  and  maximum  signal-to-average-clutter  ratio; 

1-D  spatial  frequency  (radians/unit  length); 

cartesian  components  of  2-D  spatial  frequency  Q  =  u  +  v; 

complex  conjugates  of  signal  spectra  at  the  sensor  output; 

clutter  power-spectral  densities  at  the  sensor  output; 

Fourier  transformation; 

Fourier-transform  parameters  and  convolution  variables  corresponding  to  x .  y . 


The  following  relations  with  -  denoting  the  vector  inner  product  further  define  some  of  the  quantities. 

Q2  *  Q-  Q  =  u2  +  v2  Q-  (x  +  y)  *  ux  +  vy  Q-  (a  +  0)  =  ua  +  (3a,b,c) 

Explicit  forms  of  the  power-spectral  densities  are  needed  to  evaluate  the  transforms  for  the 
impulse  responses.  The  power  spectra  employed  are 
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where  L  is  the  clutter  correlation  length.  These  equations  are  spatially  isotropic  clutter  models. 
Their  source  and  the  reasons  for  their  choice  are  discussed  in  Appendix  A.  Substituting  Eqs.  (4a.b) 
and  (3a)  into  Eqs.  (2a,b)  gives 
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The  First  terms  in  Eqs.  (5a.b)  are  the  same  as  the  impulse  responses  of  matched  filters  in  white  noise. 
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In  each  case  s  is  the  known  signal.  If  Eqs.  (6a,b)  are  differentiated  repeatedly  with  respect  to  x  and 
y,  the  even  derivatives  can  be  identified  with  the  remaining  terms  of  Eqs.  (5a,b).  The  resulting 
impulse  responses  with  inconsequential  parameters  omitted  are 
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where  V2  represents  the  Laplacian  operator  (32/3jc2)  +  (32/3y2). 
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These  analog  impulse  responses  can  be  converted  to  digital  form  by  approximating  the  deriva¬ 
tives  numerically  at  sampled  points  of  the  signals.  Appendices  B  and  C  discuss  the  approximations 
and  their  representation  by  convolution  arrays  [or),  [V2],  etc.  The  digital  matched- filter  impulse 
responses  written  in  terms  of  these  arrays  are 


[h(-0>]  -  {1  -  aL2[<J]  +  bL*[u* ]  -  clV*])  *  [s 03)]  and  (8a) 


[h(-a,  -0)]  =  11  -  aL2[V2]  +  M.4[V4]  -  cL6[ V6])  *  [s(a,  0)],  (8b) 

where  the  asterisk  denotes  convolution.  [s(/3)]  and  (s(a,  £J)J  are  arrays  of  signal  samples.  With  tabu¬ 
lated  signals,  the  impulse  responses  need  not  be  delayed  to  satisfy  causality;  hence  x  and  y  of  Eqs. 
(7a-c)  have  been  equated  to  zero.  In  keeping  with  Eqs.  (7a-c),  the  notation  indicates  that  the  opera¬ 
tors  act  on  unreversed  signals,  and  the  resulting  arrays  are  rotated  180°  about  their  centers  to  obtain 
the  impulse  responses.  The  impulse-response  arrays  are  rotated  another  180°  before  being  applied  to 
filter  input,  as  explained  in  Appendix  B.  Table  l  gives  array  approximations  of  the  derivative  opera¬ 
tors  for  equal  sampling  intervals  in  the  x  and  y  directions.  Appendix  C  explains  how  to  construct  the 
operator  arrays  when  the  x  and  y  sampling  intervals  are  not  equal. 


An  example  will  illustrate  the  use  of  Eq.  (8b)  to  determine  an  impulse  response.  To  simplify 
the  example,  it  will  be  assumed  that  the  second-derivative  term  dominates  Eq.  (8b),  allowing  the 
other  terms  to  be  ignored.  Let  the  signal  with  unit  total  strength  be 


[J<«,  0)] 


jTo  r 

6  L3  2J 


(9) 


This  could  be  a  fixed  signal  or  the  average  of  a  variable  signal.  Since  the  signal  is  known  to  be  zero 
beyond  the  listed  values,  peripheral  zeros  can  be  added  to  permit  convolution  with  the  Laplacian 
operator.  Equal  x  and  y  sampling  intervals  are  assumed  in  order  to  allow  use  of  [  V  *]  from  Table  1 . 
Then  from  Eq.  (8b)  and  Appendix  B 
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Reversing  the  resulting  array  gives  the  unnormalized  impulse  response. 


0  2  3  0 

2  -4  -10  3 

[A]  -  -  i  -2  4  0'  (11) 

0  1  0  0 

The  signal  output  from  this  unnormalized  filter  is 

0  0  2  3  0 

08  9-43 

I  6  -7  -40  -7  6 
3-4  980 

0  3  2  0  0 

with  maximum  strength  (40/6).  The  normalized  impulse  response  that  passes  the  signal  with  max¬ 
imum  strength  equal  to  the  input  total  strength  is 

0  2  3  0 

3  2  -4  -10  3 

[h]  =  “  20  1  -2  4  0 

0  1  0  0 

The  normalization,  or  gain,  of  (3/20)  is  arbitrary.  Since  gain  does  not  affect  signal-to-clutter  ratio, 
any  gain  is  permissible.  Consequently,  the  most  convenient  normalization  criterion  can  be  used  in 
each  case. 


0  2  3  ol  Ta  1 

2-4-10  3  j_  H01  l 
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0  1  0  0  [  ‘  '  °88 


For  simplicity  only,  one  term  of  the  impulse  response  was  kept  in  the  example.  In  practice  the 
terms  kept  are  determined  primarily  by  the  relative  values  of  the  power-spectral  coefficients 
(1,  aL2,  bL 4,  cL6).  These  are  best  obtained  by  fitting  the  model  Eqs.  (4a,b)  to  measured  power  spec¬ 
tra  of  the  backgrounds  involved.  Unless  one  of  the  coefficients  is  much  larger  than  the  others,  the 
impulse  response  is  a  sum  of  two  or  more  terms  properly  weighted  by  the  coefficients.  However,  the 
unity  term  is  usually  negligible  and  can  be  omitted.910  Furthermore,  unless  it  is  omitted,  the 
impulse-response  weights  do  not  sum  to  zero  and,  therefore,  the  filter  output  does  not  have  zero  mean 
locally  and  globally.  This  property  is  essential  for  correct  operation  of  the  usual  type  of  adaptive- 
threshold  sensor  that  follows  the  filter. 23-24  Omitting  the  unity  term  from  the  impulse  responses 
corresponds  in  the  frequency  domain  to  omitting  the  same  term  from  the  power-spectral  models,  and 
this  gives  transfer  functions  Eqs.  (la,b)  that  are  zero  at  zero  frequency. 

(B)  Discussion 

The  unity  term  is  sometimes  called  the  white-noise  term  because  it  has  the  same  form  as  the 
impulse  response  of  a  matched  filter  in  white  noise.  Nonetheless,  this  term  comes  from  the  constant 
of  the  clutter  model,  and  that  constant  represents,  not  white  noise,  but  power-spectral  density  of  pure 
clutter  at  zero  frequency.  The  filters  of  Eqs.  (8a. b)  are  derived  from  models  that  take  no  account  of 
white  noise.  Yet  the  clutter  waveform  actually  is  the  sum  of  a  pure-clutter  component  from  the  vari¬ 
able  background  radiance  and  a  white-noise  component  from  photon  arrival  fluctuations  and  charge 
fluctuations  in  the  detector  and  electronics.  The  true  clutter  power  spectrum  is.  therefore,  the  sum  of 
a  pure-clutter  term,  a  white-noise  term,  and  cross- spectral  terms  between  the  pure  ciutter  and  the 
white  noise.23  26  The  pure-clutter  term  is  much  larger  than  the  white-noise  and  cross-spectral  terms 
except  at  high  frequencies.  Can  the  small  terms  be  neglected?  Figure  4  of  reference  15  displays 
SNR  as  a  function  of  noise  type  for  a  mauhed  filter  in  white  noise  and  a  matched  filter  in 
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clutter.  The  SNRs  of  the  two  filters  differ  by  4.5  dB  when  both  operate  in  white  noise.  Thus, 
including  the  small  terms  in  the  (1  If2)  filter’s  design  could  not  improve  its  white-noise  performance 
more  than  4.5  dB,  and  probably  would  improve  it  much  less.  If  this  example  is  typical,  and  the  indi¬ 
cated  deficit  in  white-noise  performance  in  acceptable,  then  Eqs.  (4a, b)  are  adequate  for  deriving 
clutter-suppression  filters. 

It  is  also  worthwhile  to  consider  the  relation  between  the  clutter  models  and  the  impulse- 
response  terms  whose  coefficients  are  aL2,  bL 4,  cL6.  Comparison  of  Eqs.  (4),  (5),  and  (7)  shows 
that  or  uV  in  the  models  leads  to  jk(dk  /dyk)  or  jm  +H(dm  +n  !bxmbyn)  in  the  impulse  responses 
for  even  k ,  m ,  and  n .  Because  of  that  fact,  the  impulse-response  terms  with  coefficients 
aL2,  bL4,  cL6  are  called,  respectively,  matched  filters  for  (1  If2),  (1  If4),  and  (l//6)  clutter.  Odd 
powers  of  frequency  have  been  left  out  of  the  clutter  models  because  they  do  not  give  such  simple, 
general,  and  easily  evaluated  terms  in  the  impulse  response  (Appendix  A). 

The  isotropy  of  the  models  also  affects  the  impulse-response  terms.  Power  spectra  of  isotropic 
models  contain  the  spatial  frequencies  u,  v  only  in  the  combination  Q2  =  u2  +  v2  (ref.  27,  sec. 
10.3.2;  ref.  28,  pp.  244-248).  For  comparison,  consider  the  anisotropic  power  spectrum, 

N(u2,  v2)  =  iV0/[(l  +  L2u2)(  1  +  LyV)],  (14) 

which  leads  to  the  matched-filter  impulse  response, 


h( a,  0)  = 


-  Lf-fj  +  L2!2-^ 

ay ‘  dx-dy- 


s(x  -  a,y  -  0). 


(15) 


The  last  term  on  the  right  side  usually  predominates;  in  that  case  numerical  approximation  of 
(d*/dx2dy2)  by  [w22]  of  Table  1  gives  the  fourth-derivative  filters  of  Otazo  et  al.9-10  Equation  (15), 
unlike  Eq.  (7b),  contains  no  sixth  derivatives  and  no  homogeneous  partial  derivatives  above  the 
second.  Both  differences  are  easily  identified  with  the  powers  of  u  and  v  in  the  respective  clutter 
models.  Another  notable  difference  is  that  the  anisotropic  impulse  response  does  not  contain  the 
Laplacian  operator,  whereas  the  isotropic  response  Eq.  (7c)  can  be  written  in  powers  of  that  operator. 
This  is  entirely  due  to  the  latter  clutter  model’s  isotropy  and  the  absence  of  odd  powers  of  frequency, 
specifically  to  the  facts  that  the  correlation  length  is  a  single  constant,  and  the  spatial  frequency  com¬ 
ponents  occur  only  in  the  combination  (u2  +  v2)p ,  p  =  1,  2.  3. 


In  what  circumstances  might  one  prefer  an  isotropic  filter  or  vice  versa?  References  9-13  dis¬ 
cuss  this  and  related  questions,  so  comments  here  will  be  brief.  If  an  optical  sensor's  output  has  a 
predominant  anisotropy,  a  filter  mated  to  that  anisotropy  is  preferable.  A  sensor  could  produce  such 
an  output  if  it  views  only  a  small  class  of  anisotropic  scenes,  or  if  its  construction  or  perspective 
induce  a  strong  anisotropy  in  the  output.  On  the  other  hand,  a  sensor  may  view  a  large  class  of 
scenes  with  different  equally  probable  anisotropies,  and  it  may  be  constructed  and  operated  so  as  to 
induce  no  particular  anisotropy  in  the  output.  In  that  case  the  output  is  likely  to  be  isotropic  on  aver¬ 
age,  making  an  isotropic  filter  desirable.  Stated  somewhat  differently,  when  the  anisotropy  of  the 
sensor  output  is  unknown  and/or  highly  variable,  an  isotropic  matched  filter  could  be  a  desirable 
compromise  on  average.  However,  an  LMS  filter  may  be  preferable  in  these  circumstances  since  its 
clutter-model  parameters  change  to  suit  the  local  background. 


III.  1-D  LMS  FILTERS 


Although  2-D  filters  probably  are  used  more  often  than  1-D  filters  for  electro-optical  signal  pro¬ 
cessing,  1-D  LMS  design  is  studied  at  length  in  this  section  because  its  key  features  carry  over  to  and 
illuminate  the  more  intricate  2-D  problem,  which  is  not  easily  treated  in  the  same  detail.  Subsection 
III  A  develops  the  1-D  design  procedure  and  gives  all  that  is  needed  to  calculate  an  LMS  digital  filter 
for  any  sampled  1-D  signal.  Subsection  III  B  elaborates  certain  aspects  of  the  1-D  analysis  and 
filters. 
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(A)  Basic  Analysis  and  Example 


In  the  LMS  design  procedure  a  deterministic  model  of  signal-plus-clutter  is  fit  by  least  squares 
to  the  observed  sensor  output.  The  signal  is  represented  by  the  product  of  an  amplitude  A  and  a 
shape  function  s.  The  clutter  is  assumed  to  be  a  continuous  function  of  position  y  and  to  be 
representable  within  the  filter’s  limited  span  by  a  few  terms  of  the  function’s  Maclaurin-series  expan¬ 
sion.  With  K  samples  per  signal  length  and  a  clutter  polynomial  of  degree  M ,  the  signal-plus-clutter 
model  can  be  written 

M 

f j  —  Asj  +  £  Bmyf,  j  =  1,2,3  ...  K,  and  K  >  (M  +  2) .  (16) 

m  “0 

The  superscript  of  Bm  is  not  an  exponent;  it  simply  associates  Bm  with  a  power  of  ys  in  the  clutter 
model.  Equation  (16)  is  fit  to  the  K  observations  within  the  filter  by  adjusting  the  parameters  A  and 
Bm .  Consequently,  the  signal  amplitude  and  clutter  coefficients  vary  over  the  scene,  whereas  the  sig¬ 
nal  shape  and  orientation  are  fixed.  In  this  way  the  signal-plus-clutter  model  adjusts  to  the  input,  as 
illustrated  in  Fig.  1  which  shows  how  linear,  quadratic,  and  cubic  clutter  polynomials  might  fit  the 
background  at  three  filter  positions. 

At  any  point  within  the  filter,  the  difference  between  the  signal-plus-clutter  model  and  the 
observed  sensor  output  v;  is 

ei  mfj  ~  V  J  -  I*2*  3,  ...  K.  (17) 

Equations  (17)  are  called  equations  of  condition.  If  the  parameters  and  observations  are  equal  in 
number,  these  equations  can  be  solved  for  parameters  that  make  all  differences  zero.  Usually,  how¬ 
ever,  least-squares  equations  are  derived  by  minimizing  the  sum  of  the  squared  differences  with 
respect  to  the  parameters,  and  then  the  least-squares  equations  are  solved  for  the  best-fit  parameters. 
Either  way,  the  solution  for  the  signal  amplitude  A  provides  the  relative  weights  of  an  LMS  digital 
filter  whose  output  estimates  the  amplitude  of  the  model  signal  in  the  K  samples  contained  by  the 
filter. 


Least-squares  equations  will  now  be  derived  for  a  quintic  signal-plus-clutter  model.  To  simplify 
notation  the  subscript  k  will  imply  summation  over  the  observations,  e.g., 

ek2  *  E  Ifj  -  vy  )”•  (18' 

im  i 

Differentiating  ek  with  respect  to  the  adjustable  parameters  and  equating  the  derivatives  to  zero,  one 
obtains 


K  de,  dek  dek 

E  c.  — —  m  ek  — —  =  0  and  ek - 

/T,  ’  dA  dA  d(Bm) 


=  0,  m  =  0.  1,  2,  3.  4.  5. 


(19) 


Combining  Eqs.  (16,17,19)  gives  the  desired  least  squares  equations,  which  in  matrix  form  are 
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where,  eg.  skvk  ■  -  s;  y;  and  the  dummy  index  of  summation  on  the  right  side  has  been  changed 
from  k  to  X  for  clarity  in  later  work.  Equations  for  lower-degree  signal-plus-clutter  models  are 
obtained  simply  by  eliminating  appropriate  elements  from  Eq.  (20).  For  example,  the  equation  for  a 
quadratic  model  is  obtained  by  eliminating  the  last  three  rows  and  columns  of  the  coefficient  matrix 
and  the  last  three  elements  of  the  parameter  and  observation  vectors.  The  coefficients  of  the  parame¬ 
ters  and  observations  are  obtained  from  the  known  values  of  the  model  signal  s,  and  the  sample  co¬ 
ordinates  \j .  Both  the  coefficient  matrix  and  computation  of  its  elements  are  simplified  by  choosing 
the  co-ordinate  origin  judiciously.  A  1-D  signal  has  either  an  odd  or  even  number  of  samples  which 
lie  equally  on  each  side  of  a  central  point.  If  the  origin  is  placed  at  the  central  point,  the  number  of 
co-ordinate  values  y ,  is  halved  except  for  sign,  and  all  matrix  elements  vf  with  p  odd  are  zero.  After 
the  coefficients  have  been  determined,  the  set  of  linear  equations  is  solved  numerically  for  the  best-fit 
signal-amplitude  parameter  A  ,  and  the  terms  of  the  solution  are  written  as  a  convolution  sum  (defined 
and  illustrated  in  Appendix  B).  Finally,  the  LMS  filter  weights  are  obtained  by  inspection  of  the  con¬ 
volution  sum. 


An  LMS  filter  for  the  signal  [5  2  1  3  4  3]  determined  by  six  throws  of  a  fair  die  will  be  com¬ 
puted  as  an  example.  The  signal  is  chosen  to  emphasize  the  fact  that  the  procedure  is  applicable  to 
any  signal  whatsoever.  Assume  that  the  clutter  within  the  signal  length  is  adequately  fit  by  a  qua¬ 
dratic  function.  The  co-ordinate  origin  is  placed  at  the  signal  center,  and  v,  is  measured  in  units  of 
the  sampling  interval.  Thus  we  have: 


Model 
Indices  (j) 

Co-ordinates  (y; ) 

Signal  values  (r, ) 
Observations  (v;) 


fj  -  As,  +  (fl°  +  B 1  y,  +  Bzyf).  j 

1  2  3  4 

5>  _2  _!  1 

2  2  2  2 

5  2  13 

v,  v,  v3  v4 
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2  o 


4  3 


v5 
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Elements  are  eliminated  from  Eq.  (20)  to  suit  the  model,  and  the  co-ordinates  and  signal  values  are 
substituted  into  the  reduced  equations.  This  gives  the  following  set  of  equations  for  the  best-fit  model 
parameters: 
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The  solution  of  interest  is  that  for  the  signal  amplitude. 

A  =  (3920s xvx  -  8085vx  +  224vxvv  -  1260vv:  vx)/23856.  (22) 


The  denominator  is  the  determinant  of  the  coefficient  matrix:  it  is  absorbed  eventually  in  the  filter 
normalization  (gain).  The  unsigned  numbers  in  the  numerator  are  the  minors  of  the  elements  in  the 
first  column  of  the  determinant  array.  Since  the  denominator  is  irrelevant  and  the  subscripts  imply 
summation. 
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A  =  V  (3920s.  -  8085  4-  224v,  -  I260v.-)v.  =  V  u-,-_.,v,. 


(23) 


’  =  i 


The  right  side  of  Eq.  (23)  can  be  interpreted  as  a  convolution  sum  in  which  the  parenthetical  subsums 
are  the  LMS  filter  weights.  In  terms  of  convolution  arrays  [Appendix  B.  Eqs.  iBl-B4i)  the  last 
equality  of  Eq.  (23)  is 

A  «  [H'1w2H'3w4H'JH'6]  *  KvjVjV^jvJ  -  w6v,  +  wjVj  +  w4Vj  +  w3v4  +  w:v5  w,v#.  (24) 

By  comparison  of  Eqs.  (23)  and  (24), 

w6  -  3920*1  -  8085  +  224y,  -  I260y?  »  3080.  (25) 

In  die  same  way  one  finds  for  the  other  weights  w}  -  -3416,  w4  »  -4592.  w,  »  3472. 
w2  *  5096,  wj  “  —3640.  These  weights  can  be  checked  by  calculating  A  of  Eq.  (23)  with  the 
model-signal  values  Sj  substituted  for  the  observations  v; .  The  result  should  be  and  is  the  deter¬ 
minant  23856  of  the  coefficient  matrix.  This  check  is  based  on  the  fact  that  the  filter  output  is  pro¬ 
portional  to  the  amplitude  of  the  model-signal  content  of  the  input.  Another  check  is  based  on  the 
fact  (discussed  in  subsection  m  B)  that  an  LMS  Biter  rejects  Mac launn- senes  input  corresponding  to 
the  terms  of  its  clutter  model.  In  this  example  the  filter  should  and  does  reject  the  arbitrary  quadrat* 
input 

[(a  +  by j  +  cy*)]  «  (l/4)[(4a  -  106  +  25c)  (4a  -  66  +  9eM4a  -  26  +  c) 

(4 a  +  26  +  c)  (4a  +66  +  9c)  (4a  +  106  +  25c )[ 

Since  any  LMS  filter  blocks  constant  input,  the  weights  should  always  sum  to  zero  This  check  is 

part  of  the  previous  one,  but  is  noted  specially  because  it  is  easy  to  use.  The  filial  step  is  to  normal¬ 

ize  the  weights  according  to  any  desired  criterion,  as  discussed  in  the  matched-filter  example 

The  following  points  related  to  solution  of  Eq.  (20)  should  be  nosed. 

(1)  The  signal  values  and  observations  are  equal  in  number 

(2)  The  number  of  observations  must  equal  or  exceed  the  number  of  model  parameters 

(3)  If  the  parameters  and  observations  are  equal  in  number,  the  best-tit  parameters  can  be  calculated 
from  either  the  equations  of  condition  or  the  least-squares  equations,  otherwise  the  least -squares 
equations  must  be  used. 

(4)  With  a  nonsymmetnc  signal  the  clutter  model  may  end  with  either  an  even  or  odd  power  ot  . 
With  a  symmetric  signal  the  clutter  model  must  end  with  an  odd  power  ot  v 

The  second  and  third  points  are  familiar  mathematical  facts,  (he  fourth  receives  additional  attention  n 
the  next  subsection. 

(B)  Further  .Analysis  and  Discussion 

This  subsection  gives  further  examples,  explores  the  solution  or  Eq  20 *  in  greater  jetail  anu 
discusses  the  make-up  of  the  filter  weights  and  how  the  filters  operate  Readers  *ith  title  >r  no 
interest  in  these  matters  can  slum  or  omit  (hi>  part  Thes  ma\  *ant  :o  . onsioer  t  ater  noweser  t 
they  are  puzzled  by  certain  aspects  ot  their  calculations  tor  particular  tillers 


Numerical  solutions  like  the  preceding  one  are  easily  obtained  and  are  the  most  efficient  way  to 
derive  LMS  filters  for  given  signals.  On  the  other  hand,  algebraic  solutions  provide  more  informa- 
tkm  about  the  properties  of  the  equations  and  filters.  Accordingly,  algebraic  solutions  for  centrally 
symmetric  signals  will  be  derived  and  discussed.  These  signals  receive  special  attention  for  two  rea¬ 
son*.  First,  their  symmetry  introduces  certain  unique  features  into  the  solution  of  Eq.  (20).  Second, 
point-source  signals  are  of  this  type,  and  detection  of  infrared  point  sources  is  important  in  practice. 

With  the  co  nrrlinmr  origin  at  the  center  of  a  centrally  symmetric  signal,  all  terms  containing 
odd  powers  of  y;  are  zero  in  the  coefficient  matrix  of  Eq.  (20).  and  the  equations  separate  into  two 
seta.  One  set  contains  the  odd  clutter  parameters,  the  other  the  even  clutter  parameters  and  the  signal 

amplitude. 
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The  solution  of  Eq.  (26)  for  the  best-fit  signal  amplitude  is 


(26) 


(27) 


A  •  I J*2 1 UxVjJ  ~  I I <vx)  +  U*yt2ICyx2v,J  -  \sky£\(y£\)J.  (28) 

where  the  determinant  of  the  coefficient  matrix  is  ignored  because  it  is  absorbed  eventually  in  the 
filler  normalization.  The  symbols  |  sf 1 ,  |  st  |  .  etc.  represent  the  minors  of  the  elements  sf,  sk ,  etc. 
in  the  determinant  array,  c.g.. 


det.  of 


y*:  y*4  y? 


V4  V6  V  8 

^y*  y*  y* 


(29) 


Since  the  subscript  X  implies  summation.  Eq.  (28)  can  also  be  written 

A  *  £  l!vl*;  -  Ij*I  +  Ut.vfl.v/  -  |styt4|v/|vr  (30) 

i  *  i 

The  subsums  m  braces  are  the  impulse- response  weights  of  the  fifth-order  LMS  filter  for  the  signal 
amplitude  parameter.  Comparison  of  these  expressions  with  Eq.  (Bl)  of  Appendix  B  shows  that  the 

weights  are 

=  sk~  s.  -  !  sk  '  i  -**>'*“  I  >*,“  ~  I  Skyk  I  X*.  j  =  1.2.  ...  K.  (31) 


Least-squares  equations  applicable  with  third-  and  first-degree  clutter  models  can  be  obtained  by 
eliminating  appropriate  elements  from  Eqs  (26,27).  Correspondingly,  elimination  of  terms  from  Eq. 
(30)  gives  solutions  of  the  reduced  Eq  (26i.  Specifically.  LMS  filters  based  on  cubic  or  linear  clutter 


models  are  obtained  by  eliminating  ( |  sky?  \  y*)  or  ( |  sky{  |  yj  -  |  sky£  \  y*)  from  Eq.  (30).  The 
resulting  solutions  with  the  minors  expanded  are 

^  *  E  -  (y?)2]Sj  -  [(^)(y*4)  -  (y*2)(W)  1 

j- 1 

+  l(i*)(y*)  -  toi*y*)ly/l  V  and  (30a) 

A  -  2  (toy  “  **)  vy  (30b) 

y-i 

Again,  the  subsums  in  braces  are  the  impulse-response  weights  and  are  to  be  evaluated  as  previously 
described.  The  expanded  fifth-order  solution  is  not  given  because  of  its  length.  Specific  fifth-order 
solutions  are  best  obtained  directly  from  Eq.  (30). 

Table  2  gives  algebraic  solutions  of  Eq.  (30)  for  linear,  cubic,  and  quintic  models  and  two  sym¬ 
metric  sample  patterns.  Actually,  the  sampling  is  hardly  ever  symmetric  because  the  signal  phase  is 
random.  That  is  usually  accounted  for,  however,  by  averaging  the  sample  values  over  the  equally 
probable  phases.  With  a  symmetric  signal  the  averages  have  one  of  the  patterns  in  Table  2,  depend¬ 
ing  on  the  ratio  of  signal  length  to  sampling  interval.  On  average,  therefore,  the  examples  of  Table  2 
include  all  types  of  sampled,  1-D,  symmetric  signals. 

No  solutions  for  even-order  models  are  given  in  Table  2  because  none  exist.  The  signal  sym¬ 
metry  does  not  allow  the  clutter  function  to  end  with  an  even  power  of  .  If  one  attempts  to  use  an 
even-order  model,  he  finds  that  the  number  of  observations  must  be  sufficient  to  permit  solution  of 
Eq.  (27)  for  the  next  higher  odd  clutter  parameter.  The  fact  that  this  is  due  to  the  signal  symmetry  is 
not  apparent  in  the  least-squares  equations,  but  it  can  be  found  in  the  equations  of  condition  by 
expanding  the  coefficient  matrix  on  the  row  or  column  containing  the  signal  values. 

Table  2  provides  examples  of  several  filters  for  each  signal-plus-clutter  model.  The  first  exam¬ 
ple  for  each  model  has  the  least  number  of  observations  needed  to  determine  the  model  parameters. 
Each  successive  example  then  has  one  more  observation  than  its  predecessor.  All  filter  weights  are 
weighted  sums  of  the  model  signal’s  even  derivatives  estimated  at  the  sample  point  nearest  the  signal 
center.  This  is  illustrated  by  the  examples  of  third-order  filters,  those  based  on  a  third-degree  clutter 
polynomial. 

•  In  the  first  example,  there  are  just  enough  samples  to  estimate  the  central  value  of  the  signal's 
fourth  derivative,  and  the  filter  weights  are  multiples  of  the  estimate  {2b  -  8a  +  6)  ■  7. 

•  In  the  second  example,  the  number  of  samples  is  sufficient  to  estimate  the  signal’s  fifth 
derivative,  but  it  is  zero,  and  so  only  the  fourth  derivative  is  estimated  at  the  two  points 
nearest  the  center.  The  filter  weights  are  multiples  of  the  estimate  {b  -  3u  +  2)  «  5. 

•  In  the  third  example,  there  are  enough  samples  to  estimate  central  values  of  both  the  fourth 
and  sixth  derivatives,  and  the  filter  weights  are  weighted  sums  of  the  two.  However,  the 
derivative  values  enter  the  filter  weights  in  different  proportions  and  so  do  not  factor  out  as 
they  do  in  the  previous  examples. 

•  In  the  fourth  example,  the  filter  weights  are  weighted  sums  of  the  mode!  signal’s  fourth  and 
sixth  derivatives  estimated  at  either  point  nearest  the  signal  center. 

•  Presumably,  higher  and  higher  even  derivatives  enter  the  filter  weights  as  the  number  of  sam¬ 
ples  is  increased. 


•m 
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The  same  type  of  comments  describe  the  first-  and  fifth-order  filter  weights,  except  that  they  are 
made  from  different  derivatives. 

Figure  2  shows  that  the  derivative  components  of  the  weights  and  the  terms  of  the  clutter  func¬ 
tion  are  interchangeable  in  a  certain  sense.  For  example,  with  five  observations  and  a  linear  clutter 
function,  the  filter  weights  contain  the  second  and  fourth  derivatives.  With  the  same  number  of 
observations,  adding  quadratic  and  cubic  terms  to  the  clutter  function  eliminates  the  second  derivative 
from  the  weights.  Further  inspection  of  Fig.  2  provides  additional  examples  of  this  type.  It  follows 
that  the  weights  of  Table  2  are  built  up  from  values  of  the  model  signal’s  even  derivatives  above  the 
degree  of  the  clutter  function. 

Since  the  filter  weights  are  made  from  the  signal’s  central  derivatives,  the  weights  will  be  very 
small  or  zero  if  the  derivative  components  are  very  small  or  zero.  Obviously  such  filters  should  be 
avoided.  The  cases  in  Table  2  where  the  derivative  value  factors  out  of  the  weights  may  appear  to  be 
exceptions.  One  could  argue  that  only  the  relative  weights  matter  in  such  cases,  and  the  derivative 
Actor  can  be  ignored.  If  the  ignored  factor  is  zero,  however,  the  filter  output  is  nevertheless  zero 
with  the  signal  centered  on  the  relative  weights.  Moreover,  since  the  weights  vanish  with  their 
derivative  components,  it  can  be  seen  from  Eqs.  (31,30,28)  that  in  these  cases  the  coefficient  matrix 
of  Eq.  (26)  is  singular  and  there  is  no  filter  solution.  Thus  the  filter  always  vanishes  with  the  model- 
signal’s  central  derivatives. 

From  the  make-up  of  the  impulse  responses,  it  follows  that  the  filters  of  Table  2  operate  by 
correlating  input  with  weighted  sums  of  a  model-signal’s  even  derivatives;  however,  the  operation  can 
be  explained  more  physically  as  follows.  Let  the  input  at  any  instant  be  expanded  in  a  Maclaurin 
series  about  the  filter  center.  (This  series  represents  the  actual  signal  and  clutter;  it  is  not  related  to 
the  signal-plus-clutter  model.)  The  filter  weights  are  symmetrically  distributed  and  sum  to  zero.  As 
a  result,  the  filter  eliminates  the  constant  and  odd  terms  of  the  input  series.  Furthermore,  the  weights 
are  so  related  that  even-input  terms  below  the  last  term  of  the  clutter  model  also  are  eliminated.  Thus 
the  filter  operates  by  passing  only  even-input  terms  higher  than  the  filter  order.  This  explanation  can 
be  checked  by  working  out  examples  from  Table  2. 

The  following  useful  picture  of  an  LMS  filter’s  operation  emerges  from  the  preceding  discus¬ 
sion.  Let  the  input  consist  of  additive  signal,  pure  clutter,  and  random  noise,  each  expanded  in  a 
Maclaurin  series.  The  filter  blocks  input  terms  that  correspond  to  its  clutter  model  or  that  have  sym¬ 
metry  different  from  its  model  signal.  The  remaining  input  terms  are  correlated  with  combinations  of 
the  model-signal's  derivatives  above  the  filter  order,  and  the  result  is  passed.  The  passed  quantity  is 
a  weighted  sum  of  signal,  pure-clutter,  and  random-noise  terms  of  an  order  higher  than  the  clutter 
polynomial.  Hence,  one  can  say  that  an  LMS  filter  assigns  to  clutter  and  rejects  input  (including  sig¬ 
nal  and  random  noise)  that  can  be  represented  by  the  clutter  model  or  has  the  wrong  symmetry,  while 
it  assigns  to  signal  and  passes  input  (including  pure  clutter  and  random  noise)  that  correlates  with  the 
model-signal’s  derivatives  above  the  filter  order. 


It  is  worthwhile  to  compare  1-D  LMS  filters  for  sampled  and  continuous  symmetric  signals. 
The  least-squares  equations  and  solutions  for  the  two  filter  types  are  very  similar,  differences  being 
due  only  to  summing  or  integrating  the  squared  errors  over  the  signal  length  depending  on  whether 
the  data  are  sampled  or  continuous.  Consequently,  least-squares  equations  for  continuous  data  can  be 
obtained  from  Eqs.  (26,27)  by  substituting  signal  length  2 Y  for  number  K  of  samples  per  signal 
length,  and  replacing  sums  over  samples  by  integrals  over  signal  length.  The  same  substitutions  in 
Eqs.  (30,30a,30b)  give  the  continuous  solutions  for  the  best-fit  signal-amplitude  parameters.  For 
example,  with  a  symmetric  signal  and  a  linear  clutter  model. 

A  *  v(y)dy  -  2y{_  |.r(y )  -  <s  >  |  v(yWv  =  -v)v(v  )dy.  (32) 
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The  first  equality  is  the  continuous  analog  of  Eq.  (30b).  The  angular  brackets  <  >  in  the  second 
equality  denote  an  average  over  the  signal  length.  Inspection  of  the  last  equality  gives  the  continuous 
impulse  response, 

h(y)  =*  s(-y )  -  <s>,  -Y  £  y  £  Y 

=  0,  elsewhere.  (33) 


From  this  expression  it  is  easy  to  show  that  the  continuous  impulse  response,  like  the  discrete  one,  is 
a  weighted  sum  of  the  model  signal’s  even  derivatives  above  the  filter  order.  Replacing  the  signal 
function  with  its  Maclaurin  series  everywhere  in  Eq.  (33)  yields 
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where  (DHs) 0  is  the  model  signal’s  nth  derivative  at  the  signal  center.  The  derivative  weights  depend 
on  the  y  co-ordinate  in  the  continuous  case  and  on  the  sample  position  in  the  discrete  case.  Similar 
results  can  be  obtained  for  higher-order  filters.  All  even  derivatives  above  the  filter  order  occur  in 
continuous  impulse  responses,  but  discrete  impulse  responses  contain  no  derivatives  beyond  the 
highest  that  can  be  estimated  with  the  number  of  samples  available.  Decreasing  the  sampling  interval 
increases  the  number  of  samples  in  the  signal  length  and  brings  higher  derivatives  into  the  filter 
weights.  This  indicates  that  a  discrete  filter  passes  smoothly  to  a  continuous  one  as  the  sampling 
interval  is  reduced  to  zero. 

The  behavior  described  also  implies  that  oversampling  should  improve  LMS-filter  performance. 
Extra  samples  allow  use  of  a  higher-order  filter  or  a  more  selective  lower-order  filter  whose  weights 
contain  a  greater  number  of  derivatives.  The  lower-order  filter  passes  more  input  terms  containing 
both  signal  and  clutter,  and  is  the  better  performer  if  the  resulting  fractional  increase  of  output  signal 
exceeds  that  of  output  clutter.  In  a  staring  system  for  detection  of  point  sources,  the  needed  oversam¬ 
pling  requires  high  spatial  resolution,  several  times  higher  than  is  currently  available.  With  a  scan¬ 
ning  system  this  limitation  can  be  overcome  in  the  scan  direction  by  increasing  the  electronic  sam¬ 
pling  rate.  It  has  been  suggested  that  a  similar  effect  can  be  achieved  in  the  cross-scan  direction  by 
using  several  columns  of  detectors  with  offset  centers.29  It  is  also  possible  to  increase  the  number  of 
samples  by  adding  zeros  at  the  ends  of  the  model  signal,  but  that  is  not  as  effective  as  oversampling 
for  two  reasons:  The  zeros  provide  little  extra  information  about  the  signal,  and  the  added  clutter 
values  may  be  weakly  correlated  with  those  inside  the  nonzero  signal  length. 

Finally,  reference  15  shows  that  a  first-order  LMS  filter  and  a  matched  filter  in  (1  If1)  clutter 
have  the  same  transfer  function  for  a  raised-cosine  signal.  This  has  led  to  the  misconception  that  an 
LMS  filter  is  always  equivalent  to  some  matched  filter.  On  the  contrary,  the  two  types  of  filters  usu¬ 
ally  are  not  the  same.  Relations  that  determine  the  circumstances  in  which  they  are  identical  can  be 
obtained  by  equating  general  expressions  for  their  impulse  responses  or  transfer  functions  -  e.g..  Eqs. 
(7a)  and  (33)  for  1-D  continuous  filters  and  Eqs.  (8a)  and  (31)  for  1-D  discrete  filters.  This  problem 
is  not  studied  here  because  it  is  only  marginally  related  to  computation  of  filter  weights,  the  main 
topic  of  the  article. 

IV.  2-D  LMS  FILTERS 


The  LMS  design  procedure  can  be  applied  in  two  dimensions  almost  as  easily  as  in  one. 
although  this  has  not  been  done  heretofore.  The  derivation  of  the  least-squares  equations  is  the  same 
in  both  cases:  there  are  simply  more  equations  with  more  terms  in  the  2-D  case.  This  makes  the  2-D 


equations  more  difficult  to  solve  but  does  not  affect  the  main  features  of  the  solutions.  Since  the  1-D 
and  2-D  cases  are  so  alike,  they  are  presented  in  the  same  way.  Subsection  IV  A  contains  all  the 
information  needed  to  calculate  the  weights  of  a  2-D  LMS  digital  filter  for  any  specified  signal.  Sub¬ 
section  IV  B  provides  additional  examples  and  discussion  of  the  equations,  solutions,  and  filters  for 
symmetric  signals. 

(A)  Basic  Analysis  and  Example 


As  in  the  1-D  case,  a  model  of  signal-plus-clutter  is  fit  by  least  squares  to  sensor  output  within 
the  filter,  which  is  the  same  size  as  the  signal.  The  signal  model  is  the  product  of  an  amplitude  A 
and  a  shape  function  s,  and  the  clutter  model  is  an  arbitrary  continuous  function  represented  within 
the  filter  by  a  few  terms  of  its  truncated  Maclaurin  series.  A  right-handed  co-ordinate  system  is 
employed  with  positive  y  axis  rightward  and  positive  x  axis  downward  (Fig.  3  and  Appendix  B).  To 
simplify  notation,  the  samples  of  signal,  co-ordinates,  and  sensor  output  are  labeled  with  a  single 
ordered  subscript.  The  indexing  system  is  fully  specified  by  the  examples  in  Table  3.  With  these 
conventions  the  signal-plus-clutter  model  of  degree  M  for  K  samples  can  be  written 
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As  before,  the  superscript  of  is  not  an  exponent  but  simply  associates  the  parameter  with  the 

corresponding  powers  of  jc;  and  Vy.  A  and  B(m  _"MI  are  adjusted  to  achieve  the  least-squares  fit.  In 
operation  the  filter  algorithm  makes  these  adjustments  spontaneously  at  every  point  of  the  scene. 

To  limit  the  size  of  the  2-D  least-squares  equations,  they  will  be  derived  for  a  cubic  model 
(Af  =  3)  instead  of  the  quintic  model  used  in  the  l-D  case.  Equation  (17)  still  represents  the 
equations  of  condition  for  the  differences  e}  between  model  values  and  observations.  Differentiating 
the  sum  of  the  squared  differences  with  respect  to  the  model  parameters  and  equating  the  derivatives 
to  zero  gives 
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As  in  Eq.  (19),  the  subscript  k  implies  summation  over  the  K  samples  in  the  filter.  Equations  (17) 
and  (35)  are  now  substituted  into  Eqs.  (36)  to  obtain  the  least-squares  equations: 
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In  this  equation,  e.g.,  skxfyk  m  £  sjxfyj<  and  the  dummy  index  of  summation  has  been  changed 
from  k  to  X  on  the  right  side  to  facilitate  subsequent  analysis.  Equations  for  higher-degree  models  are 
no  more  difficult  to  derive  but  are  inconveniendy  large.  The  equations  for  quartic  and  quintic  models 
have,  respectively,  (16  x  16)  and  (22  x  22)  coefficient  matrices.  Eliminating  appropriate  elements 
from  Eq.  (37)  gives  equadons  for  lower-degree  signal-plus-clutter  models.  For  example,  eliminating 
the  last  four  rows  and  columns  of  the  coefficient  matrix  and  the  last  four  elements  of  the  parameter 
and  observation  vectors  gives  the  least-squares  equadons  for  a  quadradc  model.  Unlike  1-D  signals, 
2-D  signals  do  not  always  have  a  central  point  about  which  the  samples  are  equally  distributed.  Even 
so,  the  co-ordinate  origin  should  be  placed  as  symmetrically  as  possible  among  the  samples  in  order 
to  simplify  the  coefficient  matrix  and  computadon  of  its  elements. 

Equation  (37)  and  its  reduced  forms  are  applicable  to  any  signal  whatsoever.  They  are  solved 
for  the  filter  weights  in  the  same  way  as  the  1-D  Eq.  (20).  An  example  will  illustrate  the  procedure. 
Find  the  LMS  filter  for  the  signal  and  clutter  model  below. 


Model:  /,  =  Asj  +  [fl00  +  (fl‘%  +  fl01y;)  +  (B20xf  +  Bnx;y;  +  Bmyf)],  j  =  1,  2,  3,  ...  9 
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The  co-ordinate  values  result  from  putting  sample  1  in  the  third  quadrant  of  the  x-v  plane  (Fig.  3). 
placing  the  origin  between  the  fifth  and  sixth  samples,  and  measuring  x  and  y  in  units  of  the  sampling 
interval.  If  the  interval  is  different  in  the  two  directions,  either  interval  can  be  used  for  the  unit  of 
measure,  as  explained  below  in  subsection  IV  B.  First,  Eq.  (37)  is  reduced  to  its  quadratic  form,  and 
the  elements  of  the  coefficient  matrix  are  calculated  from  the  signal  values  and  co-ordinates 
(Xj ,  y;  ),  with  the  result 
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The  solution  for  the  best-fit  signal  amplitude  is 

A  =  (7872sxvx  -  3300vx  -  3616xxvx  -  336yxvx 


-  5232xx:vx  +  9888xvvxvx  -  14732vx: vv)/28400. 


(39) 


The  denominator  is  the  determinant  of  the  coefficient  matrix,  and  the  numbers  in  the  numerator  are, 
except  for  signs,  the  minors  of  the  elements  in  the  first  column  of  the  determinant  array.  As  in  the 
1-D  case,  the  determinant  can  be  omitted  since  it  does  not  affect  the  relative  filter  weights  and  disap¬ 
pears  finally  in  the  normalization.  Then,  taking  account  of  the  summation  implied  by  the  subscript 
and  factoring  the  common  multiple  v;  from  the  terms  of  the  numerator  gives 

9 

A  -  £  (1812s j  -  3300  -  361dc;  -  336y; 

7-1 

-  5232x/  +  988 &x;y;  -  14736y/)v;  =  £  w(10_y)vy.  (40) 

7  “  I 

The  right  side  of  Eq.  (40)  can  be  interpreted  as  a  2-D  convolution  sum  in  which  the  parenthetical  sub- 
sums  are  the  LMS  filter  weights.  This  is  most  easily  seen  if  the  last  equality  of  Eq.  (40)  is  written  in 
terms  of  convolution  arrays.  From  Appendix  B,  Eqs.  (B5-B7),  and  the  discussion  of  those  equations, 

*9  w6  0  V1  v4  V7a  0  w3  0  V1  v4  v7  a 

H>g  w5  W2  V2  V5  Vg  W{  W4  W7  V,  V5  Vg 

w7  w4  w,  v3  v6  v9  ~  w2  W  5  Wg  *  V3  v6  v9 
0  w3  0  V7  v^  0  w6  w9  Vj,,  v7  V9a 

=  W9Vi  +  WgV2  +  W7V3  +  W6V4  +  W5Vj  +  W4v6  +  W3V7  +  W2Vg  +  w,v9.  (41) 

Zero  weights  and  corresponding  unused  observations  (v^ ,  v7a,  v^)  are  included  in  the  arrays  of  Eq. 
(41)  for  consistency  with  the  notation  of  Appendix  B.  Comparing  Eqs.  (40)  and  (41)  one  finds 

w9  *  7872r,  -  3300  -  3616*:,  -  336y, 

-  5232x  ,2  +  9888*  ,y,  -  14736yf  =  -1344.  (42) 

In  the  same  way  the  other  relative  weights  are  found  to  be  =  -4384.  w7  =  5728,  >v5  *  6096. 

w5  *  -2800,  w4  =  1456,  w3  =  —4752,  w2  =  —7072,  w,  =  7072.  The  weights  are  checked  as  in 
the  1-D  case.  “<4”  calculated  from  Eq.  (40)  with  vy  =  s}  should  be  the  determinant  28400  of  the 

coefficient  matrix,  and  the  filter  should  reject  arbitrary  Maclaurin-series  input 

[(a  +  bxj  +  cyj  +  dxf  +  gxjyj  +  ry/))  corresponding  to  the  terms  of  the  clutter  polynomial.  The 

sum  of  the  weights  must  be  zero  for  rejection  of  constant  input  in  the  second  check.  The  final  step  is 

to  normalize  the  relative  weights  (set  the  filter  gain)  according  to  any  convenient  criterion. 

Everything  needed  to  calculate  an  LMS  filter  for  any  2-D  signal  has  now  been  given.  The 
reader  is  cautioned,  however,  that  a  2-D  LMS  filter  cannot  always  be  designed  simply  by  solving  Eq. 
(37).  With  quadratic  and  higher-degree  clutter  functions,  there  is  more  than  one  filter  for  a  given  sig¬ 
nal.  This  should  not  be  surprising;  a  matched  filter  for  a  given  signal  also  depends  on  the  clutter 
model,  as  Eqs.  (la,b)  show.  The  mathematical  reasons  for  the  dependence  are  different  in  the  LMS 
case,  however,  and  are  useful  for  choosing  among  alternative  LMS  filters.  This  is  discussed  in  the 
next  subsection,  where  solutions  of  Eq.  (37)  for  certain  symmetric  signals  are  considered. 


miwnmininuummiun 


(B)  Further  Analysis  and  Discussion 

The  signal  symmetries  of  interest  are  those  which  cause  all  terms  of  the  coefficient  matrix 
involving  an  odd  power  of  x  or  y  to  vanish.  Examples  of  these  symmetries  for  linear  and  higher- 
degree  signal-plus-clutter  models  are 


linear 

c  b  d 
a  1  a 
d  b  c 


higher 

degree 

c  b  c 
a  1  a 
c  b  c 


rwnn 


The  symmetry  for  a  linear  model  is  slightly  lower  because  Eq.  (37)  for  a  linear  model  contains  no 
coefficients  of  the  type  skxgyk,  p  and  r  a  1.  With  signals  of  the  specified  symmetry,  eighty-six 
matrix  elements  vanish,  and  Eq.  (37)  splits  intg  four  equations: 
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The  signal-amplitude  solution  of  Eq.  (43a),  omitting  the  determinant,  is 

*  =  I  **  1 5\vx  -  |  sk  |  vx  +  |  skxf  |  xx2  vx  -  1  skyk-  |  yx2  vx 


(44) 


■  L  1 1**2  I  sj  ~  I5*  i  +  I W  | xf  -  | skyk-  | y/)vy .  (45a) 

where  |  sf 1 ,  |  sk  |  ,  etc.  are  the  minors  of  the  elements  in  the  first  column  of  the  determinant  array: 
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Equation  (45a)  gives  the  best-fit  A  for  cubic  signal-plus-clutter  models.  As  in  the  l-D  case,  there  are 
no  solutions  for  quadratic  models  with  signals  having  the  specified  symmetry.  After  elimination  of 
the  third  and  fourth  rows  and  columns  from  Eq.  (43a),  the  signal-ampiitude  solution  (without  deter¬ 
minant)  for  a  linear  model  is  found  to  be 

K 

A  =  £  \Ksj  -  s4)  \j.  (45b) 

/  =  i 

This  has  the  same  form  as  the  corresponding  l-D  solution  Eq.  (30b). 

The  right  sides  of  Eqs.  (45a,b)  are  convolution  sums  in  which  the  filter  weights  are  the  subsums 
enclosed  by  braces.  Since  the  impulse  response  is  reversed  before  convolution  with  the  input  v; ,  the 
filter  weights  for  cubic  and  linear  signal-plus-clutter  models  are  respectively 

w(K+D-j  =  I**2 If/  ~  1**1  +  iWl*/  “  \Wk\y?  and  (47a) 

w(K  +  i)-j  ~  tej  ~  sk,  j  =  1,  2,  ...  K.  (47b) 

Expansion  of  the  minors  in  Eq.  (45a)  shows  that  all  the  terms  in  braces  are  equidimensional  in 
x ,  y ,  and  s .  That  is,  x  occurs  to  the  same  power  in  every  term,  as  also  do  y  and  s .  This  means 
that  the  relative  filter  weights  are  not  changed  if  ail  values  of  Xj ,  or  v, ,  or  are  multiplied  by  a  con¬ 
stant,  a  relation  with  three  notable  consequences:  (1)  The  signal  and  co-ordinates  can  be  scaled  to 
simplify  computations  if  the  numbers  are  awkward.  (2)  If  the  x  and  y  sampling  intervals  are 
unequal,  either  can  be  used  as  the  unit  of  measure.  (3)  A  difference  in  the  two  sampling  intervals 
affects  the  relative  filter  weights  through  the  signal  values,  not  through  the  co-ordinates.  Further 
implications  of  these  facts  are  pointed  out  in  the  last  paragraph  of  this  subsection. 

Tables  3  and  4  show  filters  for  selected  symmetric  signals  combined  with  linear  and  cubic 
clutter  functions.  A  clutter  function  is  cubic  if  it  has  a  term  xfyj  with  (p  +  r)  =  3.  The  cubic 
models  were  determined  by  eliminating  terms  from  the  complete  cubic  function  until  a  solution  for  A 
containing  all  values  of  Sj  was  found.  The  complete  linear  function  was  used  in  all  cases.  Tables  3 
and  4  begin  with  the  signals  having  the  fewest  samples  required  to  determine  the  parameters  of  a 
linear  or  cubic  model.  The  other  signals  were  chosen  to  illustrate  the  fact  that  the  signal  shape,  the 
clutter  function,  and  the  number  of  observations  jointly  determine  the  derivative  components  of  the 
filter  weights.  In  the  l-D  examples  (Table  2)  the  number  of  observations  is  varied  with  the  clutter 
function  and  signal  shape  fixed.  In  the  2-D  linear  examples  (Table  3)  the  clutter  function  is  fixed,  but 
the  signal  shape  changes  to  maintain  the  specified  symmetry  as  the  number  of  observations  is 
increased.  And  in  the  2-D  cubic  examples  (Table  4)  both  the  clutter  function  and  the  signal  shape 
change  with  the  number  of  observations.  Consequently,  the  development  of  the  filter  weights  from 
the  model-signal’s  derivatives  is  not  illustrated  as  straightforwardly  in  two  dimensions  as  in  one.  The 
key  features  are  the  same  in  both  dimensions,  however,  and  are  evident  in  Tables  2-4.  When  the 
number  of  observations  equals  the  number  of  model  parameters,  the  filter  weights  are  made  from  a 
single  model-signal  derivative.  As  the  number  of  observations  within  the  signal  is  increased,  more  or 
different  derivatives  enter  the  filter  weights  depending  on  the  signal  shape  and  clutter  function.  This 
suggests,  by  analogy  to  the  l-D  case,  that  discrete  2-D  filters  pass  smoothly  to  continuous  form  as  the 
sampling  interval  approaches  zero. 

Since  the  filter  weights  are  made  from  the  model-signal's  derivatives,  the  two  vanish  together. 
It  is  also  true,  as  can  be  seen  from  Eqs.  (47a.45a),  that  the  coefficient  matrix  in  Eq.  (43a)  is  singular 
when  the  derivative  components  of  the  filter  weights  are  zero.  Hence  the  filter  vanishes  with  the 
derivatives,  even  if  there  is  only  one  component  that  factors  out  of  the  weights,  as  in  the  first  and 
fifth  examples  of  Table  4. 


The  selectivity  of  the  filter  depends  on  the  number  of  derivative  components  in  the  weights. 
The  number  of  observations  needed  to  bring  a  given  number  of  derivatives  into  the  weights  is  usually 
greater  in  two  dimensions  than  in  one  because  of  the  greater  number  of  model  parameters  involved. 
This  makes  oversampling  of  the  signal  more  important  for  filter  performance  in  two  dimensions  than 
in  one. 

It  was  noted  earlier  that  an  LMS  filter,  like  a  matched  filter,  is  not  unique  but  depends  on  the 
clutter  model  as  well  as  the  signal.  Table  5  exhibits  alternative  third-order  LMS  filters  for  two  sig¬ 
nals  of  Table  4.  Comparison  of  the  tables  shows  that  clutter  terms  and  model-signal  derivatives  are 
interchangeable  in  the  following  way.  Adding  B^yf  to  the  second  signal-plus-clutter  model  of  Table 
4  eliminates  (D‘s  )4  -  0  from  the  filter  weights.  Similarly,  adding  B^yf  or  fl20*,2  to  the  third 
model  eliminates  (D~s )36  —  5  or  (D/s  )3  6  —  a  from  the  filter  weights.  This  effect  has  at  least  two 
consequences.  First,  in  some  instances  terms  must  be  omitted  from  the  clutter  polynomial  in  order 
for  the  filter  weights  to  contain  as  much  information  as  possible  about  the  signal.  That  is  why  the 
models  of  Table  4  were  determined  by  eliminating  terms  from  the  full  cubic  clutter  function  until  a 
signal-amplitude  solution  containing  all  signal  values  was  found.  Second,  the  filters  of  Table  5  reject 
some  input  passed  by  the  corresponding  filters  of  Table  4.  This  is  be  ise  the  filters  of  Table  5  have 
more  clutter  terms  and,  as  pointed  out  in  subsection  HI  B,  an  LMS  alter  rejects  input  representable 
by  its  clutter  function.  Which  filters  are  better  depends  on  whether  rejection  of  additional  input,  con¬ 
taining  both  signal  and  clutter,  reduces  output  signal  or  output  clutter  by  a  greater  fraction.  The 
answer  probably  differs  from  case  to  case. 

The  filter  weights  and  derivative  components  in  Tables  3-5  are  for  equal  x  and  y  sampling  inter¬ 
vals.  How  does  changing  one  of  the  intervals,  say  Ajc  ,  affect  these  results?  The  question  is  answered 
by  the  relations  noted  earlier  between  filter  weights,  co-ordinates,  and  signal  values.  With  Ay  the 
unit  of  measure,  changing  Ax  scales  Xj  values  but  does  not  alter  the  expressions  for  the  relative 
weights  as  sums  of  signal  samples  a,  b,  etc..  But  the  invariance  is  only  algebraic  because  one  or 
more  of  the  signal  values  is  changed,  thereby  altering  the  relative  weights  numerically.  At  the  same 
time,  derivative  components  involving  Ax  are  changed  by  appropriate  multiples,  and  the  coefficients 
of  these  components  in  the  filter  weights  are  changed  inversely  by  the  same  multiples,  so  that  the 
algebraic  expressions  for  the  relative  weights  as  sums  of  signal  samples  remain  unchanged.  In  practi¬ 
cal  terms  this  means  that  if  the  definitions  of  a,  /3,  etc.  are  substituted  into  the  relative  filter  weights 
in  Tables  3-5,  the  resulting  expressions  for  the  weights  in  terms  of  a.  b.  etc.  are  valid  for  any  A* 
and  Ay.  (In  fact  the  weights  were  derived  as  sums  of  signal  values  and  then,  to  show  their  basic 
structure,  were  converted  to  the  tabulated  sums  of  derivative  values  with  Ax  =  Ay  =  1  sample  inter¬ 
val.) 

V.  SUMMARY 

Spatial  filters  are  often  used  to  enhance  signal-to-clutter  ratio  in  an  electro-optic  sensor  s  output 
prior  to  signal  detection.  Matched  filters  and  LMS  filters  are  two  types  commonly  employed.  A 
matched  filter  maximizes  the  ratio  of  filtered-signal  energy  to  average  filtered-clutter  power.  An 
LMS  filter  implements  a  regression  analysis  of  the  sensor  output  and  estimates  the  amplitude  of  the 
signal  present. 

This  article  gives  methods  of  computing  the  impulse-response  weights  of  one-  and  two- 
dimensional  digital  matched  filters  and  LMS  filters  for  any  signal  and  clutter  representable  by  certain 
types  of  models.  The  methods  and  models  are  applicable  to  outputs  from  scanning  or  staring  optical 
sensors  and  to  signals  from  extended  or  point  sources.  The  signal  at  the  sensor  output  must  be  known 
to  design  either  kind  of  filter.  With  resolved  signals,  whose  shapes  are  variable,  the  design  tech¬ 
niques  are  limited  to  cases  where  an  average  signal  is  adequate  or  a  particular  form  of  the  signal  must 
be  detected.  For  matched-filter  design  the  clutter  power  spectrum  at  the  sensor  output  also  must  be 
known,  but  for  LMS  design  prior  knowledge  of  the  background  is  not  essential,  though  it  mat  be 


helpful  for  choosing  the  clutter  model.  Images  of  sensor  output  are  the  best  sources  of  the  required 
information  about  signals  and  backgrounds. 

The  matched-filter  design  technique  is  developed  by  using  isotropic  power-spectral  clutter 
models,  Eqs.  (4a,b),  which  contain  only  even  powers  of  spatial  frequency.  This  leads  to  fairly  simple 
expressions,  Eqs.  (7a-c),  for  the  1-D  and  2-D  continuous  impulse  responses  as  weighted  sums  of  the 
model  signals  and  their  derivatives.  The  weights  are  the  coefficients  of  terms  in  the  power-spectral 
models.  Consequently,  these  coefficients  must  be  known  either  theoretically  or  from  background 
measurements,  preferably  the  latter.  The  1-D  impulse-response  equation  contains  only  even  deriva¬ 
tives;  the  2-D  equation  contains  only  powers  of  the  Laplacian  operator  [(32  /  dx2)  + 
(d2  /  dy2)Y ,  p  =  1,  2,  3.  These  equations  are  converted  to  discrete  form  by  numerical  approxima¬ 
tion  of  the  derivatives.  The  results  are  a  pair  of  operator  equations  (8a,b)  for  the  filter  weights  and  a 
table  of  derivative  operator  arrays  (Table  1).  2-D  arrays  are  tabulated  only  for  equal  x  and  y  sam¬ 
pling  intervals.  The  derivations  of  the  equations  and  operators  are  given  in  sufficient  detail  so  that  a 
reader  can  repeat  them  with  different  clutter  models  and  unequal  x  and  y  sampling  intervals.  If  the 
clutter  model  contains  odd  powers  of  frequency,  the  impulse  response  is  not  easily  determined.  For 
that  reason,  such  models  were  not  used  in  this  work. 

LMS  filters  are  derived  from  least-squares  fits  of  deterministic  signal-plus-clutter  models,  Eqs. 
(16,  35)  to  sensor  output.  The  clutter  models  employed  here  are  truncated  and  sampled  Maclaurin- 
series  expansions  of  continuous  functions.  Their  coefficients,  unlike  those  of  the  power-spectral 
models,  need  not  be  known  because  an  LMS  filter  constantly  and  automatically  adjusts  the  parameters 
of  its  signal-plus-clutter  model  to  suit  the  input,  as  in  Fig.  1 .  However,  information  about  the  back¬ 
ground  is  useful  for  choosing  the  clutter  function  since  it  is  not  specified  by  LMS  design,  which 
rather  makes  best  use  of  the  model  chosen.  Least-squares  equations  (20)  are  given  for  1-D  LMS 
filters  based  on  first-  through  fifth-degree  clutter  polynomials.  Even  the  fifth-degree  equations  are 
easily  solved  with  a  hand-held  calculator  capable  of  inverting  a  (7  x  7)  matrix.  The  equations  for 

2- D  LMS  ftlters  are  much  like  those  for  1-D  filters.  There  are  simply  more  of  them,  and  each  has 
more  terms  than  its  1-D  counterpart.  Least-squares  equations  (37)  are  given  for  2-D  LMS  filters 
based  on  first-  through  third-degree  clutter  models.  Solution  of  the  third-degree  equations  requires 
inversion  of  an  (11  x  11)  matrix.  Fourth-  and  fifth-degree  equations  are  not  difficult  to  derive  or  to 
solve  numerically  with  a  computer,  but  they  are  not  given  because  they  involve  (16  x  16)  and  (22  x 
22)  coefficient  matrices. 

LMS  filters  for  specific  signals  are  most  easily  obtained  by  solving  the  least-squares  equations 
numerically,  but  the  solutions  provide  little  understanding  of  the  filters.  For  understanding,  algebraic 
solutions  are  needed.  Accordingly,  the  equations  are  solved  algebraically  for  symmetric  signals, 
which  are  of  special  interest  and  present  special  problems.  Only  two  types  of  symmetric,  1-D.  sam¬ 
pled  signals  are  possible  on  average,  and  solutions  are  given  for  both  (Table  2).  In  the  2-D  case, 
however,  the  treatment  is  necessarily  limited  to  a  few  of  the  many  possible  symmetric  signals  (Tables 

3- 5). 

Examination  of  the  algebraic  solutions  shows  that  LMS  digital-filter  weights  are  sums  of  the 
model-signal’s  derivatives  estimated  as  near  as  possible  to  the  signal  center.  Likewise,  continuous 
LMS  impulse  responses  can  be  expressed  as  sums  of  the  model-signal's  centrally  evaluated  deriva¬ 
tives,  e.g.,  Eq.  (34).  The  derivative  components  of  the  impulse  response  and  the  terms  of  the  clutter 
model  are  interchangeable;  that  is.  adding  a  term  to  the  clutter  model  eliminates  a  derivative  from  the 
impulse  response  (Fig.  2  and  Table  5).  A  continuous  LMS  impulse  response  contains  all  signal 
derivatives  not  ruled  out  by  the  clutter  model.  The  derivative  components  of  a  discrete  LMS  impulse 
response  are  farther  limited  by  the  number  of  signal  samples.  A  derivative  cannot  appear  in  the  filter 
weights  if  there  are  not  enough  samples  to  estimate  it.  As  the  sampling  interval  is  decreased,  the 
number  of  samples  in  the  signal  span  increases,  and  more  of  the  signal's  derivatives  enter  the  filter 
weights  (Fig.  2).  This  shows  how  a  discrete  LMS  filter  passes  to  continuous  form  with  decrease  of 


sampling  interval.  Improvement  of  filter  performance  by  oversampling  also  is  indicated  since  a 
filter’s  selectivity  depends  on  the  number  of  derivatives  in  its  weights. 

The  build-up  of  the  impulse  response  from  the  model-signal’s  derivatives,  except  those  ruled  out 
by  the  clutter  model,  suggests  two  simple  explanations  of  an  LMS  filter’s  operation.  Clearly  the  filter 
correlates  input  with  the  signal  derivatives,  but  a  more  physical  explanation  can  be  given  as  follows. 
An  LMS  filter  assigns  input  to  signal,  clutter,  and  residual.  Assigned  clutter  is  all  input  representable 
by  the  clutter  model.  Assigned  signal  is  the  remaining  input  that  correlates  with  the  model  signal. 
Residual  is  any  input  left  over  (none  if  the  model  parameters  and  observations  are  equal  in  number). 
Assigned  signal,  clutter,  and  residual  all  contain  actual  signal  and  clutter  in  different  proportions. 
Assigned  clutter  and  residual  are  rejected,  assigned  signal  is  passed.  Addition  of  a  term  to  the  clutter 
model  causes  additional  portions  of  actual  signal  and  clutter  to  be  assigned  to  clutter  and/or  residual 
and  rejected.  Whether  performance  improves  depends  on  whether  this  causes  actual  signal  or  actual 
clutter  in  the  output  to  decrease  by  the  greater  fraction,  a  question  that  must  be  answered  case-by¬ 
case.  Notice  that  this  physical  explanation  essentially  describes  what  a  regression  analysis  does,  in 
keeping  with  the  LMS  filter’s  mathematical  basis. 
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The  meaning  and  use  of  these  arrays  are  discussed  n  Appendices  B  and  C  They  are  not  vectors  or  matrices  except  where 
the  symbol  t  for  ‘transpose*  occurs  In  that  case  thev  are  regarded  as  row  vectors  or  as  matrices  tor  evaluating  me  evnres 
s»on  involved.  The  arrays  operate  by  convolution,  ndicated  by  an  asterisk  *  The  'Irst  step  >t  convolution  v  arrav  'ever 
sal.  Convolving  an  operator  array  with  an  array  ot  samples  estimates  a  derivative  at  he  sample  points  Operator  irri" 
must  always  fully  over\Mp  sample  arrays  This  table  fives  operator  arrays  tor  lv  -  Aa  •  sample  ntervai  .  1)  irrjs>  >r 
unequal  x  and  y  samplini  intervals  must  he  constructed  as  described  m  Appends  C 


The  first  array  is  from  forward  and  backward  difference  approximations  of  he  Jerivauve  he  second  array  •  t 
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Table  2  —  l-D  LMS  Filters  for  Symmetric  Signals 


Table  3  —  2-D,  First-order,  LMS  Filters  for  Certain  Symmetric  Signals 
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Table  4  —  2-D,  Third-order,  LMS  Filters  for  Certain  Symmetric  Signals 
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•  a.b.c.d  may  have  any  vetoes  except  thoee  which  make  the  derivative  components  vanish.  The  symmetry  of  the  signals  is  such  that  I  s.x.^y,'-  I  x.^'-O 
wMh  p.  or  r.  or  both  odd.  Theee  ere  the  conditions  t or  Eq.  (37)  to  reduce  to  Eqa.  <43a.b,c.d). 

'  ID^y'sli  •3*#~’s/  3x^3/  at  sample  poem  i  or  j.  The  denvativee  are  estimated  from  Eqs.(C2-C7)  of  Appendix  C  with  Ax  «  Ay  “  1  sample  interval  if  the 
,  the  coefficient  matrix  of  Eq.  (43a)  is  singular,  and  the  equation  cannot  ba  solved  for  a  signal- amplitude  filter. 


*  The  relative  filter  weights  ere  m  braces.  Odmeniy  the  derivative  factor  is  an  unseen  pert  of  the  filter  normalization  (gam)  constant.  The  derivative  factor  <s 
shown  here  to  emphasize  that  the  weights  are  built  up  from  the  model- signal’s  denvativea  estimated  as  near  as  possible  to  the  signal  center,  and  that  the 
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Fig.  1  —  Adjustment  of  1-D  clutter  models  to 
local  background.  This  figure  is  discussed  in  the 
first  paragraph  of  subsection  QI  A. 


Fig.  2  —  Signal-derivative  components  (D")l  ; 
of  1-D.  LMS,  digital  filters  for  symmetric  sig¬ 
nals.  This  figure  shows  the  model-signal  deriva¬ 
tives  that  make  up  the  weights  of  first-,  third-, 
and  fifth-order  filters  with  increasing  number  of 
samples/signal  length,  i.e.,  decreasing  sampling 
interval.  Comparison  of  columns  shows  that  the 
weights  are  made  of  derivatives  above  the  degree 
of  the  clutter  model  (filter  order).  Reading  up 
the  columns,  one  sees  that  more  and  higher 
derivatives  enter  the  weights  as  the  sampling 
interval  is  decreased.  Reading  across  the  rows, 
one  sees  that  adding  terms  to  the  clutter  polyno¬ 
mial  eliminates  model-signal  derivatives  from  the 
filter  weights. 
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Filter  Order  end  Degree  of  Cluttei  Function 


Fig.  3  —  Co-ordinate  system  and  numbering  of  array 
rows  and  columns,  “s"  is  model-signal  strength;  “v”  is 
observed  sensor  output  from  actual  signal  and  clutter. 
This  figure  is  discussed  in  the  last  paragraph  of  Appendix 

B. 
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APPENDIX  A 


Clutter  Models 

The  clutter  models  used  in  the  derivations  of  LMS  filters  are  deterministic,  i.e.,  represent 
unaveraged  local  clutter  properties.  These  models  are  simply  continuous  functions  that  are  fit  by  least 
squares  to  the  background  within  the  filters.  In  such  limited  regions  the  models  can  be  represented  by 
a  few  terms  of  their  Maclaurin  series  Eqs.  (16  and  35).  Maclaurin  series  are  used  for  simplicity; 
expansions  in  any  complete  set  of  orthogonal  functions  also  could  be  used.  Whatever  the  expansion, 
the  deterministic  models  conform  to  the  input;  i.e.,  the  least-squares  values  of  the  expansion  coeffi¬ 
cients  change  as  the  background  passes  through  the  LMS  filter. 

Statistical  models  are  used  for  the  matched-filter  derivations.  These  models  are  fixed;  they  do 
not  respond  to  local  background  variations;  rather,  their  parameters  are  adjusted  before  use  to  approx¬ 
imate  certain  globally  averaged  properties  of  the  anticipated  clutter.  The  remainder  of  this  appendix 
discusses  three  candidate  statistical  models  and  the  reasons  for  choosing  one  of  them. 


Temporal  stationarity  and  spatial  isotropy  are  important  qualities  related  to  the  choice.  A  pro¬ 
cess  is  stationary  and  isotropic  if  its  average  properties  depend  only  on  time  and  distance  between  its 
points  (ref.  27,  sec.  10.3.2).  Most  types  of  optical  clutter  are  naturally  nonstationary  and  anisotropic. 
Moreover,  these  qualities  are  affected  by  experimental  conditions.  Nonstationarity  is  increased  if  the 
sensor  views  different  backgrounds  in  rapid  succession.  Anisotropy  may  be  induced  or  altered  by  the 
sensor’s  construction,  and  is  always  affected  by  the  sensor’s  perspective  because  the  length  scale 
varies  with  direction  in  the  background  plane  unless  the  scene  is  viewed  along  a  normal.11  Neverthe¬ 
less,  clutter  often  is  represented  by  stationary  isotropic  models  for  reasons  of  convenience  and 
ignorance:  Nonstationarity  and  anisotropy  are  mathematically  difficult  to  treat,  and  too  little  is  known 
about  them  to  make  the  effort  worthwhile.  In  addition,  unless  a  certain  anisotropy  predominates 
because  of  the  sensor's  construction,  clutter  will  tend  to  be  isotropic  on  average  if  the  sensor  views 
many  scenes  from  all  perspectives.12,13 

Clutter  usually  is  modeled  in  the  time  and  space  domain  by  autocorrelation  functions,  and  in  the 
frequency  domain  by  their  Fourier  transforms,  power  spectra.  In  this  article  time  dependence  and 
hence  stationarity  are  ignored.  Three  purely  spatial  autocorrelation  functions  that  have  been  proposed 
as  models  are9-13 

/f(Ar)  =■  <r  exp  ( -  |  Ar  |  /£.).  (Ah 

/?(Ax.  yAv)  =  <r  exp  [-(Aar  +  y:Av:)1/2/L].  and  (A2) 


/?(Ajc  .  Ay )  =*  <r  exp(  -(  1  Sx  /Lx )  -  ( 1  Ay  • /Lx  )J.  (A3i 

Here  Ar  =  Ax  •+■  Av  is  a  displacement  vector.  <r  is  clutter  variance.  L  is  correlation  length,  and 
is  a  factor  that  accounts  for  a  change  of  scale  due  to  perspective.  The  first  model  is  isotropic;  the 
second  is  anisotropic  due  to  perspective  alone:  the  third  is  intrinsically  anisotropic,  i.e..  without 
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regard  to  perspective  and  even  if  Lx  =  Ly.  The  power  spectra  derived  from  these  autocorrelation 
functions  are11 

=  F[tf(Ar)]  =  2mrLz/(\  +  I2!)2)3'2,  (A4) 

S{u~.  v2/^2)  =*  F[R(Ax,  7Ay)]  =  2t<t2(Z.2/7)/[1  +  L2u2  +  L2(v/7)2]3/2,  and  (A5) 

N(uz,  v2)  =  F[/?(Ar,  A y)]  =  +  IXV)(  1  +  l/v2)],  (A6) 

where  F  indicates  Fourier  transformation  and  0  =  5  +  v  is  2-D  spatial  frequency  with  x  and  y  com¬ 
ponents  u,  9.  Note  that  Eq.  (A4)  is  a  function  of  the  combination  Q2  =  (u2  +  v2),  not  a  function  of 
the  separate  variables  u2,  v2.  This  results  from  the  isotropy  and  consequent  circular  symmetry  of  the 
autocorrelation  function,  Eq.  (Al)  (ref.  11,  app.  B;  ref.  28,  pp.  244-248).  A  1-D  power  spectral 
model  is  obtained  heuristically  from  Eq.  (A4)  by  substituting  u  =  0  and  replacing  v  by  u  for  nota- 
tional  clarity. 

For  reasons  already  explained,  the  isotropic  power  spectral  model  is  used  here.  However,  Eq. 
(A4)  and  its  1-D  counterpart  are  mathematically  awkward  for  deriving  matched  filters.  More  suitable 
forms  are  obtained  by  expanding  the  denominators  in  Maclaurin  series.  This  leads  to 

(V(u2)  =  iV0/(l  +  aL2ur  +  bL*w*  +  cL 6gj6)  and  (A7) 

NiQ2)  =  N0/(  1  +  aL2&  +  bLW  +  cL^W),  (A8) 

where  N0  is  power-spectral  density  at  zero  frequency.  Terms  beyond  the  fourth  are  dropped  because 
they  lead  to  filters  larger  than  any  expected  correlation  lengths.  (N0,  aL2,  bL 4,  cL6)  can  be  regarded 
as  disposable  parameters  to  be  determined  from  measurements  of  natural  backgrounds.  In  reference 
30  an  approximation  to  a  1-D  power  spectrum  is  derived,  rather  than  heuristically  inferred,  from  Eq. 
(A4);  the  approximate  function  varies  as  1/(1  +  dw2).  Additional  terms  in  the  denominator  of  Eq. 
(A7)  can  be  regarded  as  empirical  higher-order  approximations.  Equation  (A8)  is  an  isotropic  model 

because  it  depends  only  on  Q  =  +vM2  +  v2  and  therefore  its  transform,  the  autocorrelation  function, 
depends  only  on  the  magnitude  of  A r  (ref.  28,  pp.  244-248). 

An  isotropic  power  spectrum  like  Eq.  (A8),  but  containing  both  odd  and  even  powers  of  fr,  can 
be  obtained  by  expanding  the  denominator  of  N(Q)  =  N0/f(Q)  in  a  Maclaurin  series: 

N(Q)  =  N0/{1  +  <5fl  +  efi2  +  ffi3  + . )  (A9) 

This  model  can  decrease  less  rapidly  than  (1/fi2)  as  some  measured  power  spectra  do.  Since  physical 
quantities  are  real,  however,  autocorrelation  functions  and  hence  power  spectra  must  be  real  and 
even.  This  ordinarily  is  considered  to  rule  out  odd  powers  of  frequency  in  series  or  polynomial 
power-spectral  models  (ref.  31,  sec.  5.10).  But  if  such  models  are  functions  of  frequency  magnitude, 
they  can  contain  odd  powers  and  yet  be  even.  Thus  models  with  odd  powers  of  frequency  seem 
physically  permissible.  When  Eq.  (A9)  is  substituted  into  Eq.  (2b),  however,  the  odd  powers  lead  to 
terms  in  the  impulse  response  which  cannot  be  related  simply  and  generally  to  the  signal.  This  disad¬ 
vantage  of  the  odd  powers  is  considered  to  outweigh  the  greater  flexibility  they  provide  for  fining 
measurements.  For  that  reason  models  with  odd  powers  of  frequency  are  not  used  in  this  article. 

Attention  has  been  directed  to  certain  consistencies  and  discrepancies  between  physical  proper¬ 
ties  of  clutter  and  mathematical  properties  of  statistical  models.  Though  such  comparisons  are  signifi¬ 
cant  and  useful,  it  should  be  remembered  that  the  Wiener-Khintchine  type  of  models  employed  here 


arc  primarily  mathematical  constructs  and  tools.  They  were  not  discovered  by  analysis  of  physical 
processes,  and  their  connection  with  natural  phenomena  often  is  obscure.32  Thus  their  suitability  is  to 
be  judged  more  by  results  obtained  with  them  than  by  physical  consonance  or  dissonance  between 
model  and  phenomenon. 
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APPENDIX  B 

Convolution  Arrays  and  Co-ordinates 

The  operation  of  nonrecursive  digital  filters  can  be  represented  by  tables  of  sampled  input  data, 
impulse-response  weights,  and  output  values.  In  Appendix  C  it  is  shown  that  the  same  type  of  tables 
can  be  used  to  represent  numerical  approximations  of  derivative  operators.  For  reasons  that  will 
become  apparent,  these  tables  are  called  convolution  arrays.  Reference  28  (pp.  30-40)  explains  the 
basis  of  this  representation  and  an  alternative  using  vectors  and  matrices.  The  vector  and  matrix 
representation  is  discussed  at  length  in  references  8,33,34.  Convolution  arrays  are  not  vectors  or 
matrices  because  they  do  not  multiply  as  vectors  and  matrices  do;  the  multiplication  rule  is  described 
below.  However,  addition  and  multiplication  by  a  scalar  are  performed  as  with  vectors  and  matrices. 
To  distinguish  the  two  types  of  quantities,  convolution  arrays  will  be  enclosed  in  brackets  [  ],  vectors 
and  matrices  in  parentheses.  Convolution  arrays  lack  the  mathematical  completeness  and  power  of 
vectors  and  matrices  but  are  much  simpler.  In  this  article,  convolution  arrays  are  used  wherever  pos¬ 
sible. 

The  elements  of  a  1-D  convolution  array  are  written  in  a  row  and  labeled  with  indices  j .  For 
example, 

[v,  v2  v3  . . .  \j  vj  +  l  \j+2...}  and  [w,  w2  w3], 

with  \j  and  wj  being  input  values  and  impulse- response  weights  respectively.  An  input  and  impulse 
response  are  convolved  by  rotating  the  weight  array  180°  about  its  center,  passing  it  stepwise  over  the 
input  array,  and  at  each  step  summing  the  products  of  coincident  elements.  Graphically,  a  single  step 
of  the  process  is 


[w3  W2  w,]  - 

x  x  x  (Bl) 

2 

[v,  v2  v3  . . .  \j  vy  +  l  \j+2  ...]  =  £  "3-jVj^. 

J-0 


Rotating  the  weight  array  180°  reverses  one  of  the  sequences  as  convolution  requires.  Either  array 
could  be  rotated,  but  if  the  input  sequence  is  reversed,  some  thought  may  be  necessary  to  obtain  the 
output  in  the  right  order.  Reversing  the  weight  array  avoids  the  need  for  thought.  A  quantity  like 
that  on  the  right  side  of  Eq.  (Bl),  in  which  one  index  decreases  while  the  other  increases,  is  called  a 
convolution  sum  since  it  approximates  a  single  value  of  a  convolution  integral.  Passing  the  weight 
array  over  the  input  gives  a  point-by-point  numeric  estimate  of  the  convolution  integral’s  functional 
form.  The  resulting  filter  output  is  a  row  of  convolution  sums  (ref.  28.  p.  31).  With  »  denoting 
convolution,  the  output  obtained  by  passing  a  five-element  input  through  a  two-weight  filter  is 

(w,  w;]*[v,  v2  v3  v4  v5]  =  (B2) 

[(1V2V|  T  VV,V-,)(VV2V2  +  H',V3)(W,V3  +  VV|V4)(W:V4  -I-  U'  |  v?)|. 
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The  number  n  of  input  samples  must  equal  or  exceed  the  number  m  of  weights,  and  the  weights  must 
always  fully  overlap  the  input.  Consequently,  the  output  has  (m  -  1)  fewer  values  than  the  input, 
i.e.,  n  —  (m  —  1)  values.  Sometimes  the  input  can  be  extended  by  adding  peripheral  zeros,  but  only 
when  it  is  known  to  be  zero  beyond  the  listed  samples.  If  the  input  above  were  known  to  begin  and 
end  with  zeros,  then  the  filter  output  would  be 


[W]  wj]  *  [0  v,  v2  v3  v4  v5  0]  = 

[(H',V,)(M'2V1  +  WlV2)(w2V2  +  M'1V3)(H'2V3  +  W',V4)(H>2V4  +  W ,  Vj)(w2Vj)]. 


If  two  arrays  are  the  same  size,  they  can  be  multiplied  by  placing  one  over  the  other  without  rotating 
either,  and  then  summing  products  of  coincident  elements.  Multiplication  is  indicated  simply  by  writ¬ 
ing  the  arrays  in  sequence.  Thus 


lVl  v2  V3][h»j  w2  w3]  =  W(Vj  4-  w2v2  +  W3V3. 


The  elements  of  a  two-dimensional  convolution  array  are  written  in  square  or  rectangular  form 
and  can  be  indexed  in  various  ways.  A  common  method  is  to  label  them  with  row  and  column 
numbers  ij ,  for  example. 


vn  vl2 
V21  v22 


wll  w  12 
W2,  w  22 


This  input  and  impulse  response  are  convolved  by  rotating  the  weight  array  180°  about  its  center, 
passing  it  stepwise  horizontally  and  vertically  over  the  input  array,  and  at  each  step  summing  products 
of  coincident  elements.  Graphically,  a  single  step  of  the  process  is 


*22  '*'21  *  * 

'*’12  “'ll  X  X 


Vf  1 1 


'1,+  lly  vC,-r|M/flr 


I  I 

^  1  2— a  1 1  2 —  li  i—a  *  —  ti  I  • 

a— 0  li—0 


Here  also,  rotating  the  weight  array  180°  about  its  center  reverses  it  as  convolution  requires.  The 
quantity  on  the  right  side  of  Eq.  (B5)  is  a  2-D  convolution  sum,  and  the  filter  output  is  an  array  of 
such  sums.  The  output  from  a  four-element  filter  with  a  nine-element  input  is 


1  v" 

vi: 

v13  1 

Wll 

Wp 

1 

-  *  Vi, 

Vii 

Vii  I 

w:i 

w22 

— ’  1 

V31 

v?: 

v?3 

(WiiVu  4-  W;|V)2  4-  W|iV,|  -I-  w,,v22)  (vv:2v,;  4-  w2|V,3  4-  WpVii  4-  u-nv,3) 
(WnVij  4"  WijVn  4"  IVpV-jj  4“  W]jV"p)  ( VV’ tVn  4-  W 1 1 V  i  1“  U’pVp  — 


As  in  the  1-D  case,  the  two  arrays  being  convolved  must  always  fully  overlap.  Thus,  if  the  impulse 
response  has  m  rows  and  n  columns,  the  output  array  has  (m  -  1)  fewer  rows  and  (n  -  1)  fewer 
columns  than  the  input  array.  Peripheral  zeros  may  be  added  to  extend  2-D  input  arrays  only  if  they 
are  known  to  be  zero  beyond  the  listed  samples.  The  procedure  for  multiplying  2-D  arrays  of  equal 
size  is  analogous  to  the  procedure  for  1-D  arrays. 

Elements  of  2-D  arrays  also  can  be  labeled  with  a  single  ordered  subscript.  Such  a  system  for 
matrices  is  described  in  references  8,33,34.  In  that  system,  labeling  begins  at  the  top  left  element  of 
the  array,  goes  down  the  column,  then  to  the  top  of  the  next  column,  and  continues  thus  to  the  bottom 
right  element.  Labeling  of  course  does  not  affect  the  rules  for  addition,  multiplication  by  a  scalar, 
multiplication,  or  convolution.  For  example, 

n't  w4  w7  vt  v4  v7  9 

w2  wj  wg  v2  v5  v8  =  £  WjV,-.  (B7) 

w3  w6  w9  |_V3  Vs  V9  '  =  ‘ 

This  is  not  the  only  possible  system  of  labeling  with  a  single  ordered  subscript.  The  purpose  of  such 
a  system  is  to  simplify  a  discussion  or  derivation  by  simplifying  notation,  or  to  extend  1-D  mathemat¬ 
ical  techniques  to  two  dimensions.  Any  system  that  does  either  is  permissible. 

Input  data  in  a  convolution  array  are  related  as  well  to  a  co-ordinate  system.  Consequendy,  the 
connection  between  the  arrangement  of  the  data  in  the  array  and  in  the  co-ordinate  system  should  be 
specified.  The  3-D  co-ordinate  system  used  here  is  shown  in  Fig.  3.  It  is  a  right-handed  system  with 
the  positive  x  and  y  directions  downward  and  rightward  respectively.  This  choice  of  co-ordinates 
reconciles  two  conventional  notations:  x  corresponds  to  row  number,  y  to  column  number;  x  pre¬ 
cedes  y,  and  row  number  precedes  column  number;  x  and  row  number  increase  downward,  y  and 
column  number  increase  rightward.  To  simplify  computations,  the  co-ordinate  origin  is  placed  as 
symmetrically  as  possible  among  the  samples  of  signal.  With  the  usual  rectangular  sampling  grid,  the 
origin  is  located  at  the  center  of  the  smallest  rectangle  enclosing  the  signal,  and  the  unit  of  measure  is 
the  sampling  interval  in  either  the  x  or  y  direction.  If  the  two  intervals  are  the  same,  as  in  Fig.  3, 
the  co-ordinates  of  the  sample  points  are  integers  or  odd  half  integers.  Row  and  column  numbering 
begins  with  11  at  the  third-quadrant  comer  of  the  smallest  rectangle  enclosing  the  signal.  Alterna¬ 
tively,  the  co-ordinates,  signals,  and  observations  can  be  labeled  with  a  single  subscript  in  any  con¬ 
venient  order,  as  already  explained. 


APPENDIX  C 


Approximation  of  Derivatives  and  Derivative  Operators 

This  appendix  describes  estimation  of  derivatives  by  finite  differences  and  approximation  of 
derivative  operators  by  convolution  arrays.  The  finite-difference  formulae  are  needed  to  identify 
derivative  components  of  LMS-filter  weights  and  to  derive  the  operator  arrays,  which  are  needed  to 
calculate  matched-filter  weights  from  Eqs.  (8a, b). 

1.  Finite-Difference  Approximations 

The  basic  relation  between  finite  differences  and  derivatives  is  (ref.  35,  art.  18) 

=  f«)(y  +  dnAy),  0  <  9  <  1.  (Cl) 

Ayn 

This  relation  says  that  if  an  interval  is  divided  into  n  equal  subintervals  of  length  Av ,  the  n  th  differ¬ 
ence  of  a  continuous  function  f(y)  in  the  interval,  divided  by  A y",  is  equal  to  the  nth  derivative  of 
the  function  somewhere  in  the  interval.  Equation  (Cl)  is  not  completely  satisfactory  for  applications 
because  it  does  not  specify  the  point  where  the  derivative  is  evaluated.  However,  most  numerical- 
analysis  texts  derive  more  specific  relations  from  interpolation  formulae,  e.g.,  the  Newton-Sterling 
central-difference  formula,  and  the  Gregory-Newton  forward-  and  backward-difference  formulae.35-43 
These  formulae  yield  the  following  lowest-order  approximations  of  the  first  six  derivatives  (ref.  35, 
art.  48;  ref.  36,  secs.  35,36). 

A y(dv/dy)  =  Av,  =  Vv,  =  -v,  +  v2 
Av2(d2v/dy2)  =  A:v,  =  V:v3  =  52v:  =  v,  -  2v,  +  v3 

A y2(d2\  /d\2)  =  AJv,  =  V3v4  =  — v,  +  3v,  —  3v3  +  v4  (C2-C7) 

Ay4(d4v/d'y4)  =  A4v,  =  V4v5  =  64v3  =  v,  -  4v,  +  6v3  -  4v4  +  v5 

Ay5(d5\/dy5)  =  A5v,  =  V5v6  =  -v,  +  5v;  -  10v3  +  I0v4  -  5v5  +  vh 

A v6(<f6v/dy6)  =  A6v,  =  v6v7  =  S6v4  =  v,  -  6v,  +  15v3  -  20v4  +  15v5  -  6vs  +  v- 

Here  v  is  position  or  time,  v  is  sensor  output,  Vj  is  a  sample  of  sensor  output,  and  An .  V  n .  &n  are 
forward,  backward,  and  central  differences.  The  forward-,  backward-,  and  central-difference  formu¬ 
lae  are  the  same  for  even  derivatives,  but  for  odd  derivatives  the  lowest-order  central-difference 
approximations  are 

A  y(d\/dv)  = /i6v:  =  +  v3). 

Ay2(di\ /dv2)  =  ^oiv,  =  v,  +  2v:  -  2v4  i-  v<).  (C8-C10) 

Ay5id5v/dy5)  =  *i65v4  =  ~<-vi  +  Jv:  ~  5v,  -  5v<  -  4vh  -  v-i. 
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If  the  unit  of  y  is  one  sample  interval.  Av"  =  1  and  the  formulae  as  written  approximate  the  deriva¬ 
tives  themselves.  The  central-difference  formulae  are  usually  more  accurate  than  the  forward-  and 
backward-difference  formulae  (ref.  36.  sec.  29).  The  formulae  for  the  fifth  and  sixth  derivatives  are 
included  for  completeness;  ordinarily  they  are  not  used  because  differences  higher  than  the  fourth  tend 
to  be  erratic  (ref.  35,  art.  17). 

As  noted,  these  formulae  are  the  lowest-order  approximations  of  the  derivatives.  Higher-order 
approximations  are  available  but  are  seldom  used  for  filters,  possibly  because  additional  samples  are 
required,  and  reduced  correlation  between  the  added  and  original  samples  may  offset  the  greater  accu¬ 
racy  of  the  higher-order  formulae.  In  addition,  the  higher-order  approximations  require  higher-order 
differences,  which  usually  are  avoided  as  already  explained.  Second-order  central-difference  approxi¬ 
mations  of  the  first  six  derivatives  are 

A y(d\/dy)  =  ^Vj  -  ~  8v,  +  8v4  -  v,>. 

A y~(d}\/dy2)  =  <52v3  — +  16v:  ~  30v,  +  16v4  -  v,). 

A V3(d3v/dv3)  s  #i53v4  — ^-j*£5v4  =  — (V,  -  8v2  +  I3v,  -  I3v,  +  8v6  -  v-). 

4  8 

A y\d*v/dy4)  =  5\  -  — 5*v4  =  — (-v,  +  12v;  -  39v,  +  56v4  -  39v,  -  12v6  -  v.).  (Cl  1-C16) 
6  6 

A y5(d5v/dy5)  *  ii53v3  — ^-#i57v5  *  -j(v,  -  9v:  +  26v,  -  29v4  +  29v6  -  26v-  +  9v,  -  v,). 

3  6 

Av^d^v/rfy6)  =  56v,  -  — 5*v,  =  — ( -v,  +  I2v,  -  52v,  +•  1 16v4  -  I50v,  1 16v6  -  52v-  +■  I2v„  - 

4  4 

Extensive  tables  containing  coefficients  of  v;  in  formulae  like  Eqs.  (C2-C16)  can  be  found  in  refer¬ 
ence  44.  These  tables  provide  many  additional  approximations  with  estimates  of  the  errors  involved 
in  their  use. 

Equations  (C2-C16)  approximate  derivatives  at  different  sample  points  from  sets  of  values  begin¬ 
ning  with  the  same  sample  v,.  The  forward-,  backward-,  and  central-difference  formulae  estimate  the 
derivatives,  respectively,  at  the  first  (v , last  (v:  to  v7).  and  central  iv:  v3.  v4,  v5»  sample  points. 
This  must  be  taken  into  account  when  approximating  partial  derivatives  and  sums  of  derivatives 
There  are  numerous  alternatives  to  these  formulae  as  indicated  by  the  tables  of  reference  44  That 
should  not  be  surprising  since  this  is  a  problem  in  approximation,  and  such  problems  usually  can  be 
solved  in  many  ways.  In  particular,  the  derivatives  need  not  be  estimated  only  at  first,  last,  and  cen¬ 
tral  points  or  even  at  sample  points:  the  interpolation  formulae  permit  their  estimation  at  arbitrary 
points.  However.  Eqs.  (C2-C7)  are  most  satisfactory  for  the  purposes  of  this  article.  If  other  formu¬ 
lae  are  used,  it  is  sometimes  impossible  to  identify  the  derivative  components  of  the  LMS  filter 
weights.  The  probable  cause  of  this  is  pointed  out  in  the  next  part. 


2.  Derivative-Operator  Arrays 

Arrays  approximating  continuous  derivative  operators  will  now  be  obtained  by  writing  Eqs. 
(C2-C10)  in  terms  of  convolution  arrays.  It  is  important  to  remember  that  the  finite-difference  equa¬ 
tions  determine  how  and  where  the  derivatives  are  estimated.  The  arrays  do  not  affect  the  approxima¬ 
tions;  they  are  simply  a  convenient  way  of  writing  the  equations. 

With  Ay  *  1  sample  interval,  the  lowest-order  central-difference  formula  for  the  first  derivative 
is 

dv/dy  =  ~[—  1  0  l]tv,  v,  v3]  *  -j(  -  v,  +  v3). 

The  row  of  weights  approximates  the  operator  (d  !dy),  but  it  is  not  very  useful  because  it  can  operate 
only  once  on  an  array  of  its  own  size.  Its  utility  is  much  increased  by  reversing  the  weights  to  form  a 
convolution  array  and  writing 

(d/dy)  =  -j[l  0  -11. 

This  array,  when  convolved  with  any  other  of  equal  or  greater  size,  operates  repeatedly  on  the  other 
array  to  estimate  the  derivative  at  the  sample  points  that  coincide  with  the  array's  central  weight,  i.e.. 

y[l  0  -1]  *  [...  ...  )  =[...  (dv/dy  )„  (dv/dv),^  .  . .  ]. 

Thus  both  the  derivative's  values  and  its  functional  form  are  approximated  numerically.  Convolution 
arrays  approximating  the  other  monovariate-derivative  operators  can  be  constructed  in  the  same  way. 
The  results  are  given  in  Table  1  where  the  array  for  the  nth-derivative  operator  is  represented  by 
[w"].  Since  these  operators  are  actually  zero  beyond  the  listed  weights,  they  can  be  extended  by 
peripheral  zeros  when  necessary. 

Before  taking  up  partial  derivatives,  it  is  necessary  to  consider  some  properties  of  the  1-D 
operator  arrays  and  implications  of  those  properties.  It  may  seem  that  the  result  of  convolving 
discrete  operators  should  be  predictable  by  analogy  to  continuous  operators.  That  is  true  for  discrete 
operators  derived  from  Eqs.  (C2-C7V- the  forward-,  backward-,  and  even  central-difference  formulae 
For  example,  convolving  [w1]  and  (w3)  from  Eqs.  (C2)  and  <C3)  gives  (w3)  from  Eq.  <C4».  as 
expected  by  analogy  to  (d/dy)  (d2/dy:)  =  (di/dyi). 

[w‘l  *  [w:]  =  [1  -  11  .  (0  I  -2  1  0|  =  [1  -3  3  -1)  =  (w3]. 

But  operators  derived  from  Eqs.  (C8-C16)— the  odd  central -difference  and  higher-order  formulae— do 
not  behave  likewise.  This  emphasizes  that  discrete  operators  represent  finite-difference  formulae  and 
must  be  used  consistently  with  them,  not  by  analogy  to  continuous  operators.  Conversely,  the 
discrete-operator  behavior  indicates  that  Eqs.  (C2-C7)  are  analytically  related  m  the  same  way  as  con¬ 
tinuous  derivatives,  but  Eqs.  (C8-C16)  are  not.  This  explains  why  the  former  equations  are  more 
satisfactory  than  the  latter  for  identifying  the  derivative  components  of  LMS-filter  weights.  In  appli¬ 
cations  where  numeric  accuracy  is  more  important  than  analytic  properties.  Eqs  (C8-Clb»  should  be 
better. 

Turning  now  to  bivariate  partial  derivatives,  the  associated  operators  are  approximated  b\  square 
or  rectangular  arrays  whose  horizontal  and  vertical  directions  correspond,  respectively,  to  ihe  \  and  \ 
variables  (Fig.  3).  The  first  step  in  constructing  a  bivariate  operator  is  to  derive  an  approximation  for 
the  partial  derivative  from  ftnue-difference  equations.  The  monovariate  parent  equations  Jctermtne 
the  point  where  the  partial  derivative  is  estimated.  Monovariate  operators  are  then  combined  to  form 


a  2-D  array  that  acts  on  data  to  reproduce  the  approximate  formula.43  The  procedure  is  best 
explained  with  examples.  Initially  it  will  be  assumed  that  Ax  *  Ay  *  1  sample  interval.  Fre¬ 
quently.  however,  the  sampling  interval  is  different  in  the  two  directions,  and  the  difference  must  be 
taken  into  account  as  described  at  the  end  of  this  appendix. 

Central-difference  operators  for  even  homogeneous  derivatives  are  easily  constructed.  For 
example,  from  Eq.  (C3) 

{dh/dy2)  a  v,,  -  2vi2  +  vx3  and  (dh/dx2)  a  v,y  -  2v2y  +  v3y,  (C17) 

with  the  derivatives  estimated  at  the  central  points.  Now  suppose  we  have  the  data  array 


Vl! 

vt2 

v13 

V21 

v22 

v23 

V31 

v32 

v33 

and  we  want  to  estimate  (A^/Ay2)  and  (A^v/Ax2)  at  the  position  of  vr.  Add  peripheral  zeros  to  [w:J 
of  Table  1  and  to  its  transpose  [w2}',  and  write  tentatively 


(A2/Ay2)  a 

0  0  o' 

l  -2  1 

-  [w02],  (A2/Ax2)  3 

0  i  o’ 
0-2  0 

0  0  0. 

.0  1  0 

Next,  apply  these  operators  to  the  data  array,  and  compare  the  results  with  Eqs.  (07). 


0  0  0 

VU  V12  V13 

1  -2  1 

♦ 

V2I  v22  v23 

0  0  0 

V31  v32  v33 

© ' 

o 

vll  VI2  v13 

0-2  0 

* 

v:i  v22  v23 

0  10 

j_v3l  v32  v33 

v:i  ~  2VJ2  +  vy  a  (A^v/Ay2)^ 


v |2  -  2v:,  +  v32  a  (A^v/Ax2)^ 


The  operators  give  the  desired  derivative  estimates.  Since  the  estimates  are  made  at  the  same  point 
and  the  arrays  are  equidimensional .  the  operators  can  be  added  to  obtain  the  approximate  Laplacian 
operator. 


0  1  0 

V2  «  (A^/Av)  +  (A2 /Ax2)  a  (w02]  +  [w20]  *  !  I  -4  1  =  [  V -]. 

i  i 

10  1  0  _ 

It  should  be  noted  that  the  peripheral  zeros  are  not  essential  in  the  arrays  for  iA2/Av2)  and  iA:/A.t:). 
which  also  can  be  written 

1  ’ 

<A*'Av>  s  (1  -2  1|  and  iA*  Sx~)  s  -2 

I 


However,  these  operators  cannot  be  added  to  obtain  the  Laplacian  operator  since  they  are  not  equidi¬ 
mensional  and  do  not  necessarily  estimate  the  derivatives  at  ihe  same  point  Convolving  an\  ot  ihese 


operators  with  a  data  array  of  equal  or  greater  size  estimates  a  derivative  or  sum  of  derivatives  at 
sample  points  that  coincide  with  the  operator  center.  For  example. 


0  1  0 
1  -4  1  • 
0  1  0 


The  question  marks  indicate  that  the  Laplacian  derivative  is  indeterminate  at  the  sample  points  in  the 
first  and  last  rows  and  columns  because  the  operator,  when  centered  at  these  points,  lies  partly  outside 
the  data  array. 

Forward-  and  backward-difference  operators  are  constructed  in  the  same  way  as  central- 
difference  operators.  For  example,  from  Eq.  (C2) 

(<N/dy)  at  -\xi  +  y,2  and  (dv/dx)  a  -vly  +  v2v,  (C18) 

with  the  derivatives  estimated  at  either  the  first  or  last  sample  points.  From  [w1]  and  [w1]'  of  Table 
1  we  have  tentatively 

(d/dy)  a  [1  -1]  -  [w01]  and  ( d/dx )  a  =  [w10]. 

To  determine  if  and  where  these  operators  approximate  the  derivatives,  apply  them  to  a  (2  x  2)  data 
array,  and  compare  the  results  with  Eqs.  (C18). 

v„  v12  (-vu  +  v,2)  (dv/dy)i,  (3v/3v),2 

^  *  v21  v22  -  ( — v2j  +  \22)  =  (d\/dy)2l  ~  (dv/3y)22 

=  t(-vu  +  v21)(-v12  +  v22)]  a  [Ov/dx),,  (3v/3x),2)] 

s  [(3v/3x)21  (3v/3x)22] 

The  operators  do  indeed  estimate  the  derivatives  at  the  indicated  points.  The  positions  of  the  esti¬ 
mates  are  ambiguous,  but  the  ambiguity  is  in  the  forward-  and  backward-difference  formulae.  The 
arrays  correctly  represent  those  formulae. 

Operators  for  mixed  partial  derivatives  are  constructed  by  the  same  methods  used  for  homogene¬ 
ous  derivatives.  Let  us  find  the  forward-  and  backward-difference  operator  for  ( d\/dxdy ).  From 
Eqs.  (C2.C3) 

(d2v/dy2)  a  vtl  -  2vt,  +  vl3, 

(d2v/dxdy~)  =  (d/dx  )(3‘v/3y)  5  (—  vn  -r  v2[)  —  2(  —  vl2  +  \22)  +  (— v,3  +  v:3).  (C19t 


? . ?\ 


VI1  • 

•  V1 J 

i  (V2v)2i . 

(VV)!0.„ 

v,l  • 

•  v</ 

i  (v2v)(<_i>2.. 

..(V  Jv)((_i)o-n 

with  the  third  derivative  estimated  at  either  vu  or  v23.  The  data  array  that  contains  only  the  samples 
in  Eq.  (C19)  is 


V11  v12  v13 


V21  v22  v23 


An  array  of  numbers  is  needed  which,  when  reversed  and  multiplied  with  this  data  array,  gives  the 
right  side  of  Eq.  (C19).  By  inspection  the  required  array  is 


[w12]  = 


1  -2  1 
-1  2  -1 


If  the  1-D  forward-  and  backward-difference  operators  [w1],  [w2]  are  regarded  temporarily  as  row 
vectors,  the  array  above  is  the  matrix  product  [vv'j'fw2].  Thus  the  desired  2-D  forward-  and 
backward-difference  operator  is 


(d3/dxdy2)  s  [w1]'^2]  =  [w12]. 


The  derivative  is  estimated  at  either  the  top-left  or  lower-right  comer  of  the  array,  the  ambiguity 
again  originating  in  the  forward-  and  backward-difference  formulae. 


Next  let  us  find  the  central-difference  operator  for  the  same  derivative.  From  Eqs.  (C3.C8) 

02v/3y2)  s  Vjj  -  2v,2  +  vx3. 


(d/dxK&v/dy2)  &  y(-vn  +  v31)  -  (-v12  +  v32)  +  y(-vl3  +  v33). 


(C20) 


with  the  third  derivative  estimated  at  the  position  of  v22.  The  right  side  of  Eq.  (C20)  in  array  nota¬ 
tion  is 


1 

-2 

1 " 

1 - 

>~ 

rt 

>" 

> 

1 _ 

0 

0 

0 

* 

V21  V22  Vij 

.-1 

2 

-l . 

- » 

ro 

> 

n 

eo 

> 

> 

_ 1 

so  the  first  array  must  be  an  operator  which  approximates  (d3\/dxdy2)  at  the  array  center.  This  array 
is  the  matrix  product  of  the  1-D  central-difference  operators  [w1]'  and  [w2]  regarded  temporarily  as 
row  vectors.  Thus  the  desired  2-D  central-difference  operator  is 


0 3/dxdy2)  s  [w‘]'[w:]  =  [w12]. 


This  is  the  same  as  the  relation  for  the  corresponding  2-D  forward-  and  backward-difference  operator. 
In  general,  for  any  mixed  partial  derivative. 

0"*"  /dxmdyn)  a  KnlV',l  =  [vvm,,l. 
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a 


4 


,v 


where  [w"],  [w"]  are  l-D  operators  regarded  temporarily  as  row  vectors,  and  t  denotes  the  tran¬ 
spose.  The  l-D  operators  may  represent  either  central-difference  formulae  or  forward-  and 
backward -difference  formulae.  The  formulae  determine  the  points  where  (w'’*”]  estimates  the  deriva¬ 
tive. 


Thus  far  it  has  been  assumed  that  Ay  *  Ax  =  1  sample  interval.  If  the  sample  interval  is  not 
the  same  in  the  two  directions,  the  finite-difference  formulae  must  be  expressed  in  a  common  unit 
before  the  operator  arrays  are  constructed.  Suppose  the  Laplacian -operator  array  with  At  =  3Ay  is 
required.  Either  sample  interval  can  be  chosen  as  the  unit,  but  to  avoid  fractions  in  the  arrays,  let  At 
be  the  unit.  Then  from  Eq.  (C3) 

•^•(d^v/dy2)  a  vl(  -  2vtJ  +  v3  and  (dlv/dx1)  a  vlv  -  2v:v  -+-  v3v. 


The  corresponding  operator  arrays  are 


(d*/dy2)  a 


0  0  0 
9  -18  9 
0  0  0  ! 


0  10 

(w02!  and  {dr/dxz)  a  0  -2  0  «  |w:ol. 

[o  1  0  _ 


Since  these  arrays  are  equidimensionai.  have  the  same  umts.  and  estimate  the  derivatives  at  the  same 
point,  they  can  be  added  to  obtain  the  Lap  1  ac lan-ope rator  array  for  sampled  data  with  At  =  3Av 


[V21  =  (w02]  4-  [w'°] 


ro  i  o" 

9  -20  9 
0  1  0 


It  should  be  noted  that  use  of  a  common  unit  for  the  r  and  v  variables  does  not  resolve  all 
issues  raised  by  unequal  x  and  v  sampling  intervals.  In  those  circumstances,  for  example,  the 
Nyquist  frequency  and  the  accuracy  of  the  derivative  estimates  are  different  in  the  two  directions 
What  do  these  facts  imply?  Such  issues  are  worth  considering,  but  they  do  not  affect  calculation  of 
filter  weights  and.  therefore,  are  not  addressed  here. 
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